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introduced molecular-level chemistry models that predict equilibrium and 
nonequilibrium reaction rates using only kinetic theory and fundamental molecular 
properties are extended in the current work to include electronic energy level 
transitions and reactions involving charged particles. These extensions are shown to 
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Chapter 1 : Introduction 

The aerothermodynamic environments experienced by hypersonic vehicles are 
characterized by strong shock waves, which produce high temperatures resulting in 
high heating of the vehicle. The severity of the environment is dependent upon the 
mission profile and the vehicle configuration[l]. A high-speed planetary-entry blunt- 
body vehicle, such as the Galileo probe entering Jupiter or a high speed Earth return 
capsule, encounters an environment dominated by high enthalpy effects, which 
include elevated levels of ionization and may result in significant radiative heating. 
The aerodynamic and aerothennal loads on a hypersonic vehicle set the performance 
requirements for various sub-systems, such as the thermal protection system (TPS) 
and the control system. The design and safe operation of these vehicles, therefore, 
require adequate definition of the flight environment. 

During the design phase, the aerothermodynamic performance of a hypersonic 
vehicle is defined by use of either experimental or computational tools. However, as 
in most flight regimes, these definitions cannot be solely obtained from ground test 
facilities, as no facility can reproduce all aspects of the flight environment. This 
limitation is particularly true for the hypersonic flight regime because of the high 
energy requirement needed to create a true hypersonic environment at a reasonable 
scale on the ground, and especially applicable to the transitional/rarefied portion of 
the trajectory because of the very low density requirement. Therefore, designers rely 
on computational tools for the prediction of the aerothennal loads experienced by 
vehicles entering or re-entering the atmosphere from space. 

In rarefied flows, because of the breakdown of the continuum assumption, kinetic 
methods are employed. Since its introduction in the early 1960’s, the direct 
simulation Monte Carlo (DSMC) method of Bird[2] has become the de facto standard 
for simulating low-density gas dynamics. For these types of flows, traditional 
computational fluid dynamics (CFD) methods are invalid because assumptions made 
in developing the differential equations on which they are based break down under 
rarefied conditions. In contrast to continuum CFD methods, the DSMC method 
performs a direct physical simulation of the gas at the molecular level that closely 
mimics the Boltzmann equations. In the simulation, particles are tracked in space and 
time, accounting for both gas surface interactions and intermolecular collisions. 

In the upcoming years, the ability to compute rarefied, ionized hypersonic flows will 
become increasingly important. Missions such as Earth (re)entry, landing high mass 
payloads on Mars, and the exploration of the outer planets and their satellites will 
continue to push the state of the art in DSMC simulations. These missions may result 
in shock-layer temperatures reaching into the tens-of-thousands of degrees, where the 
treatment of electronic energy levels and ionization processes becomes increasingly 
important. The proper treatment of these mechanisms is needed to accurately predict 
overall shock wave thickness and structure, as well as heat transfer to the vehicle. 
Even when the reactants are minor species, the overall flow field is not affected by 
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details of the energy modes and reaction model, but a particular radiation signature 
may be strongly affected[3]. Either way, the accurate treatment and prediction of 
electronic energy distribution and ionization is an important component of the design 
of vehicle sub-systems. 

1.1 Objectives of this Work 

Although there has been much recent interest in the modeling of ionization 
processes [4-9], the inclusion of electronic energy levels in DSMC simulations has 
largely been overlooked and they are not included in many DSMC implementations. 
However, recent advances in DSMC algorithms[10-14], specifically the Quantum- 
Kinetic (Q-K) Chemistry Model, have allowed the prediction of equilibrium and non- 
equilibrium reaction rates using only kinetic theory and fundamental molecular 
properties (i.e., no macroscopic reaction-rate information). Models such as these are 
important especially when reaction rate data are limited or the gas is far from 
equilibrium and an engineering approximation is required. Thus far, these new 
models have only covered neutral-species reactions. It is the purpose of this thesis to 
expand the Q-K model to include electronic energy level transitions and also to 
include reactions involving charged species. 

In addition to ionization, radiation can become an important factor depending on the 
mission and thermal protection system. For example, inflatable ballutes can be used 
to decelerate spacecraft while still high in the atmosphere. The material used for the 
skin of the ballute would have a very low tolerance for heating, necessitating the 
knowledge of what the radiative heating load would be. Several methods of 
calculating radiative heating for rarefield flows[15-18] have been developed. For the 
current study, however, an uncoupled methodology used in many computational fluid 
dynamics radiation evaluations[19, 20] was adopted. While this study does not 
improve on the methodologies used for calculating the radiative heat transfer, having 
separate species non-equilibrium electronic temperatures would make radiative 
predictions much more accurate. 

The models developed herein will be validated against reaction rates reported in the 
literature and applied to the second flight of the Flight Investigation Reentry 
Environment (FIRE II) Project. This was an Apollo-era flight experiment to measure 
the radiative and convective heating on spacecraft material during atmospheric 
reentry at lunar return speeds. The trajectory point at 76.42 km was chosen for the 
computations because, while the flow is nearly in the transitional regime and DSMC 
is more appropriate, computational fluid dynamics (CFD) has frequently been applied 
at this altitude with success. Therefore, the DSMC simulations presented in this 
thesis will be compared to CFD solutions and flight data from this trajectory point. 
However, while the methods developed herein are completely applicable to higher 
altitudes, the use of CFD is not recommended because of the break down of the 
assumption of continuum flow. 



1.2 Outline of Thesis 


The remainder of this thesis is organized as follows. The second chapter provides an 
overview of the governing equations that are the basis for the simulation of fluid 
flows from the continuum through rarefied regimes. Also, a description of the 
methodology and structure of the DSMC algorithm, both “traditional” DSMC and the 
more recent algorithmic advances, is presented as well as previously existing 
electronic energy level treatments and ionized gas models. Chapter Three describes 
the extensions to the Q-K model and their implementation. In Chapter Four, these 
extensions are applied to an adiabatic heat bath and sampled transition/reaction rates 
are compared to those rates found in the literature. In Chapter Five, these extensions 
are then compared against previous chemical reaction models as applied to the second 
Flight Investigation of Reentry Environment (FIRE II) experimental]. Finally, 
Chapter Six summarizes the results of the thesis, the contributions of the thsis to the 
state-of-the-art, and recommended future research directions. 



Chapter 2: Simulation of Rarefied and Weakly Ionized Reentry Flows 

In this chapter, the theoretical basis of particle simulation is first summarized. Next, 
a description of the direct simulation Monte Carlo numerical method is given, 
including the traditional methods and more recent algorithmic advances, as well as 
previous attempts to include electronic energy levels and models of ionized flow. 
Throughout this description, mention will be given of any modifications that were 
required to be performed on the DSMC algorithm of choice. 

2.1 Theoretical Basis of Particle Simulation 

The particle simulation method has been developed independently by plasma 
physicists and aerothennodynamisists. Plasma physicists, who simulate charged 
particles, have developed the particle-in-cell (PIC) method based on the theory of 
charge assignment and force interpolation. This theory has made it possible to obtain 
smooth distributions of charge density and current density in Maxwell’s equation of 
electromagnetic fields. The details of the PIC method can be found in the classical 
works of Birdsall and Langdon[22] and Hockney and Eastwood[23]. On the other 
hand, aerothermodynamicists have been developing a particle simulation method of 
neutral species beginning with the work of Bird[2, 24], This method is called the 
direct simulation Monte Carlo (DSMC) method. 

Particle simulation of both charged and neutral species is necessary to examine the 
structures of high-speed, rarefied, reacting flows. However, for flows where the 
plasma density is considered low ( — 10%), the assumption of charge neutrality is 
generally assumed and the charged particles are treated as neutral particles[25-30]. 
More will be discussed on these techniques in Section 2.2.4. The current section 
describes how the probability laws used in particle simulation may be derived[3 1] 
from a relevant Boltzmann equation for a weakly ionized plasma and simplifications 
will be pointed out for our assumption of charge neutrality. 

2.1.1 The Boltzmann Equation 

The theoretical basis of particle simulation can be explained using a weakly ionized 
plasma consisting of electrons (e), ions (A), and molecules ( B ). If Coulomb collisions 
e-e, e-A, and A -A can be disregarded, only the remaining e-B, A-B, and B-B collisions 
need to be considered. The particle motions of species e, A, and B are governed by 
the respective Boltzmann equation for each of these species. Except for Coulomb 
collisions, the Boltzmann equation for electrons or ions is linear, but the Boltzmann 
equation for molecules is nonlinear because of the B-B collisions. Although the PIC 
method was developed for charged species and the DSMC method for neutral species, 
the idea behind the two methods is similar. The Boltzmann equation of ions, A, is 
written as: 
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Here, Coulomb collisions are neglected and only A-B collisions are considered, F j is 

the force acting on particle A with a mass of niA,fA{ v ) and /«( w ) are the velocity 
distribution functions of species A and B, v' and w' are the post-collision velocities, 
g - |v - w| is the relative speed, <J U1 is the differential cross section, and dQ (=sin/ 

d% d i//) is the solid angle. In the remainder of the derivation, f A (v,x,t) and 
f B (v,x,t) are expressed as f A (v) and /„(v) in the right hand side of the equation 
for simplicity. The differential cross section <J 4B is a function of g and the deflection 
angle %. The force F is given by: 


F = q A (E + vxB) (2.1.2) 

where qA is the charge of the particle A and E and B are the electric and magnetic 
fields, respectively. However, for our purposes, the force experienced by the charged 
particles is effectively zero since there is an assumption of no charge separation (i.e., 
charge neutrality). For the sake of completeness, though, this tenn will be kept in the 
remainder of the derivation. 

If the method to obtain f A (v,x,t ) is given, the solution off a can be found at any time 
by step-by-step calculations started from an initial velocity distribution function. For 
a small At, we have: 


f A (v,x,t + At) = f A (v,x,t) + At 


5A 

dt 


(2.1.3) 


Substituting df A / dt of Equation (2. 1 . 1) into Equation (2.1 .3) gives us: 

f A (v,x,t + At) = (l + AtJ){\ + AtD) f A (v,x,t) (2.1.4) 


where D and J are operators defined by 


D(f„) 


d/i ! d_ 

dx,- m. dv. 
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(2.1.5) 


j (Ja ) = JJ [/a i v ')f B i w ') - f A i v )f B ( w )] set AB dQ. dw (2.1 .6) 

D and J denote non-collisional and collisional operators, respectively, and therefore, 
Equation (2.1.4) decouples the non-collisional and collisional parts of the motion. 
Equation (2.1.4) is equal to Equation (2.1.3) if the term of order (At) 2 is disregarded 
such that the solution obtained from Equation (2.1.4) is first-order accurate. In 
general, At should be much less than the mean free time x (the mean free path divided 
by the mean speed of the particles). This condition should be strictly satisfied in the 
case when particle A is charged [32]. Note, however, that if particle A is a molecule, 
the condition can be replaced by a weaker one of At < x. The reason why the weaker 



condition can be applied is that molecular trajectories are straight and, therefore, 
accurate for any At and that the state of the gas changes little in the mean free time. 


Equation (2.1.4) suggests that a method of solving the Boltzmann equation is to first 
obtain f A = (1 + A tD)f A (non-collisional) and then obtain (l + A tj) f* A (collisional). 
The first procedure is equivalent to solving: 
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(2.1.7) 


by use of the initial distribution function f A (v,x,t ) . If the solution at t+At is 
f* A (v,x,t ) , then the solution is f* A (y,x,t ) = f A [v - [Pi m A ^At,x - vA , and 
therefore 


f A (v + Av,x + A x,t) = f A (v,x,t) 


(2.1.8) 


where A v = (E / m A ) At and Av = vAt . Equation (2.1.8) implies that the particles of 
species A move according to the equation of motion 


m A^; = q A( E + vxB ) 

dx _ _ 
dt 


(2-1-9) 


These equations of motion are used in PIC simulations which means that the 
Boltzmann equation without the collision term is that solved by the PIC method. 
When Equations (2.1.9) are solved for all particles, which is trivial for DSMC since 

there is no net force, we have f A (v,x,t) . The second procedure is equivalent to 
solving 


^ = (2.1.10) 

using the initial distribution function f A (v,x,t) . The solution of Equation (2.1.10) at 
t+At is f A [v,x,t + At) . 

Equation (2.1.10), which does not include the term df A / dx ] , represents the 

collisional relaxation of the velocity distribution in a spatially uniform state of species 
A. Equation (2.1.10) can be applied only to a small cell where the gradient of density, 
temperature, and flow velocity of species A can be regarded as uniform. Because of 
this requirement, the cell dimension should be of the order of the mean free path of 
species A, and is the basis for the cell size of the DSMC method. 



A solution of equation (2.1.10) can be obtained by an equivalent stochastic 
process[33]. First, in particle simulation, it is sufficient to consider the motions of a 
set of particles randomly sampled from all particles. The sampled particles are 
referred to as simulated particles. The number of simulated particles in a cell is much 
smaller than is the number of real particles. The weight, W, is defined as the number 
of real particles represented by each simulated particle. After moving all simulated 
particles by At, we identify the simulated particles in each cell by examining the 
positions of the particles. Let Na be the number of simulated particles in a cell and 
v A1 ,v A2 be the particle velocities. The velocity distribution function f* A (v,x,t) for 
a cell located at x is expressed as 


n a 

fl{v,x,t) = ^j-J J 5 i (v-v Ai ) (2.1.11) 

^ A i = 1 

where n A is the number density of the cell, S 3 ( ) is the Dirac delta function, and v Ai is 
the solution of Eq. (2.1.9) at time t+At. The function S' (v - v Ai ) is the probability 
density function for the velocity v Ai of particle A,. Similarly, the velocity distribution 
function of species B at the end of the particle’s motion is 

fl{w,x,t) = ^^S i (w-v B \ (2.1.12) 

7V B j~ 1 

where v Bi is the velocity of simulated particle Bj. The solution of f A (v,x,t + At) can 
be obtained from 

f A (v,x,t + At) = (l + A tJ) f* (v,x,t) 

(2.1.13) 

= fl (v,x,t) + At JJ[/ A (v')f B ( w)-f A (v)f B (w)]g(J^ dQ dw 

Substitution of Equations (2.1.1 1) and (2.1.12) into Equation (2.1.13) yields 

n a 

f A {v,x,t + At) = -^-^f Ai (v) (2.1.14) 

where f Aj (v) is the probability density function for the velocity of particle A, at t+At. 

It is given by 


A (v) = (l ■ - P^S 1 (v - v A ) + P A fi Ai (v) (2.1.15) 

The structure of Eq. (2.1.15) states that P A i is the collisional probability of particle A t . 
Note that if P Ai = 0, f Ai (v) = S ' (v - v Ai ) , and there is no change in velocity because 
of a collision. The probability P Ai is given by 



where 


P\,Bj ~ N B n B8Ai,Bj®T 


(2.1.17) 


Here, g AiRj = | v Ai - v Rj | is the relative speed and o' H is the total collision cross 

section. Equation (2.1.16) reveals that Paubj is the probability of particle A, colliding 
with particle Bj in At. The function Q Aj (v) is the probability density function for the 
post-collision velocity v . It is written as: 


e,,(v)= 



(2.1.18) 


Here, the ratio PAi.BjIPAi represents the conditional probability that the collision 
partner of particle A, is particle Bj under the condition that particle A, collides in At. 
Here, Qalbj (f ) is the probability density function for the post-collision velocity of 
particle A t when its collision partner is particle Bj, which is given by: 



(2.1.19) 


where 
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( 2 . 1 . 20 ) 


M a = m t /(mA+in/j), M B = mB/(mA+m B ), X is the angle between g' and v Ai - v Bj , and 
<3 AB (gAi,Bj>%) is the differential cross section. From Eq. (2.1.19), we have 


Qa, <5 (l g I - g Ai Bj ) dg 


AB 
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(2.1.21) 


where dv = M 3 B dg ' = M \g r dg' sin %d %dy/ . Note that Ml is the Jacobian of the 
transfonnation from v to g ' , and then Eq. (2.1.20) yields the post-collision velocity 
Ta, (= f ) of particle A, 
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(2.1.22) 
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where g'(= v' M - v' Bj ) is the post-collision relative velocity and |g'| = g Ai Bj . The 
probability that g' is included in the solid angle sin jr/jr/i// is 

® {§Ai,Bj^x) s ^ n X^X^¥ / °V {sAi.Bj) . The direction {X’V) of 8 can be sampled 
using this probability. 

Equations (2.1.16) and (2.1.17), which are derived from the Boltzmann equation, can 
also be obtained by use of the elementary free-path theory. The collision frequency 
of particle A t is 


v 


A, 



where gAi.B is the relative speed between A, and some particle B, and the bar denotes 
the average of particles B. Therefore, the average of all simulated particles B in a cell 
and the collision probability of particle A t take the fonn 


= Et A t = 


nAt -h A 


N 




AB 


B j= 1 



This equation agrees with Equation (2. 1.16). We cannot obtain the probability laws 
to detennine the post-collision velocities from the free-path theory, however. In 
many cases, we can write such probability laws intuitively because the laws derived 
from the Boltzmann equation agree with the physical image of the collision. 

The probability laws used to determine the velocity of particle A , at t+At are then: 

1) Particle Aj collides in At with a probability of Pau 

2) Under the condition that particle Aj collides, its collision partner is particle Bj 
with a probability of Pai.bJPai- 

3) If the collision partner is Bj, the post-collision velocity v' Ai is given by 
Equation (2.1.22). 

Therefore, the particle simulation method can be derived from the Boltzmann 
equations of ions and is true for any species. Although the Boltzmann equation for a 
simple gas is nonlinear because of collisions between like molecules, the simulation 
method was derived in 1 980[34] . Before this work, those in the field of rarefied gas 
dynamics had regarded the particle simulation method as a kind of numerical 
experiment. By use of the measure theory, mathematicians [3 5] verified that Nanbu’s 
particle simulation method is an appropriate method for solving the Boltzmann 
equation. 


o 



The most important ideas of the particle simulation method can be summarized as: 1) 
decoupling of collisionless motion and collisions; and 2) calculating collisions 
independently in each cell. These basic concepts have been shown to naturally result 
from the operator splitting in the form of (\+AtJ)(\+AtD). Now consider the 
conditions for choosing At. The operator splitting is allowed if and only if 
At\j(f A )| <sc f A and At |d(/ a )| <sc f A . The first condition states that the effect of 
collisions in At must be small; therefore the collision probability P Ai - At / r is small, 
where x is the mean free time. The time step must satisfy the condition At <SC T . The 
second condition requires that |Ax| <tc |x| and |Av| <s: |v| in Eq. (2.1.8). That is, 

v At <sc x and (f / m A ) At <sc v , where v is the mean particle speed and F / m A is the 

mean acceleration. The condition that v At <sc x is usually replaced by the Courant 
conditions At < A c , where A c denotes the cell size. Therefore, a particle’s 

displacement in At must be less than the cell size. The condition [f / m A j At <sc v 

requires that the change in the particle’s speed in At must be much less than the 
original speed for PIC simulations. Therefore, a smaller At should be used for an 
external field. 

2.2 Traditional Direct Simulation Monte Carlo Methods 

The DSMC technique was developed primarily by Graeme Bird, now Emeritus 
Professor at the University of Sydney, Australia, nearly a half century ago. Since its 
inception, this technique has been used in many different applications ranging from 
the simulation of micro-electromechanical devices (MEMs), to various industrial and 
pharmacological deposition simulations, to the simulation of (re)entry and other 
rarefied flow fields. 

In the work outlined herein, the DSMC method is used to simulate the Boltzmann 
equation for rarefied/transitional reentry flows in which the Navier-Stokes equations 
are not valid throughout the entire flow field. Specifically, the algorithm employed is 
the DSMC Analysis Code (DAC)[36, 37], which was developed by a team led by 
LeBeau beginning at the NASA Langley Research Center and continuing at the 
NASA Johnson Space Center. The DSMC algorithms in the release version of DAC 
are based on those developed by Bird[2], 

Regardless of the DSMC code used, the basic structure of the algorithm is the same. 
The following is included to briefly describe the flow of a generic DSMC algorithm. 

2.2.1 Collision Algorithm 

In each computational cell, particles are chosen to collide in such a way as to 
reproduce the required macroscopic collision rate. The simplest implementation of 
collision pair selection is to disregard particle position within the cell when selecting 
the pair. This methodology restricts the size of a computational cell in that it must be 
smaller than the local average mean free path of the flow field. This was the first 
selection methodology and is still utilized today. The next selection methodology 



was to select the first particle randomly in the cell and then select the particle closest 
to the first excluding its previous collision partner, or the so-called virtual sub-cell 
(VSC) algorithm of LeBeau[38]. This methodology allows for a larger collision cell 
size as long as the resulting mean collision distance divided by mean free path 
remains at least as small as the resulting value from the original selection 
methodology. Other selection methodologies have been proposed[39, 40], but the 
two discussed above are the most prevalent. The VSC algorithm is that which was 
used in the current study. 

In many modem DSMC codes, Bird’s No Time Counter (NTC) method is used to 
select colliding pairs[2]. At each time step and in each computational cell, the 
number of particle pairs to be tested for collisions is calculated as follows: 

N, M =i„N(o T g)^At (2.2.1) 

where n is the number density of gas particles in the cell, N is the average number of 
simulator particles in the cell, At is the simulation time step, g is the relative velocity 
of a colliding pair of particles, and Ot is the total collision cross section of the 
colliding pair. Each of these pairs are then selected to collide if the probability: 

p = (2 2 2) 

coll ( _ \ 

\ t / max 

is less than a random number uniformly distributed between zero and one. This 
method assumes that the collision rates between species will be implicitly reproduced, 
even though all species are lumped into the same group for the pair selection process. 
This method can become inefficient if the product of the collision cross section and 
relative velocity becomes very much different between species pairs, as is the case 
when electrons are introduced into the flow, because the maximum value of a T g is 
very large and results in many null collisions. This is the basic algorithm used by the 
DAC code. 

In order to address this issue, the collision algorithm in DAC was modified as 
follows[31]. The collision algorithm can be divided into two groups. One includes 
collisions between similar particles, and the other collisions between dissimilar 
particles. The collision term of the Boltzmann equation is nonlinear for like 
collisions and linear for unlike collisions. Consequently, the collision algorithm for 
like collisions is different from that for the unlike collisions. 

There are five classes of collisions involving charged particles: electrons with 
neutrals, ions with neutrals, electrons with ions, ions with ions, and electrons with 
electrons. The interactions of the last three classes obey the Coulomb force law: 

F _ 1 Mi 
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(2.2.3) 



where e 0 is the permittivity of free space and r is the distance between the particles. 

Bird[41] and Taylor et. al. [25], however, used the variable hard sphere (VHS) model 
of Bird to describe the first four classes of collisions by setting the diameter of the 
heavy ions equal to the diameter of their neutral counterparts. The reference diameter 
of the electrons is set equal to le-10 m, which is orders of magnitude larger than the 
classical value of 6e-15 m. This choice was somewhat arbitrary, but Bird did mention 
that the electron diameter was varied in his simulations by a factor of three with no 
significant effect on the flow field parameters. Neither modeled electron-electron 
collisions. This is the approach that will be followed herein. 

In general, then, one can divide the collision algorithm into two for weakly ionized 
plasmas, neglecting Coulomb collisions: 1) collisions between unlike particles; and 
2) collisions between like particles. The algorithm for type (1) is derived from the 
linear Boltzmann equation as in Section 2.1.1, and that for type (2) is from the 
nonlinear Boltzmann equation. 

Let us consider the elastic collisions: 

A{v a ) + b {v b )-> A{V a ) + b {v' b ) 

If gas B is at rest and in equilibrium at temperature Tb, the probability that particle A 
elastically collides with some particle B is: 


P con = n Bg° T {g) At (2-2.4) 


where ns is the number density of gas B, g = \v A - v /; | is the relative velocity, and the 
average is taken over the Maxwellian distribution of v B . The average in 
Equation(2.2.4) is for all simulated particles in the cell, which gives: 


P con = ^J J SPr(s j ) A t 


(2.2.5) 
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where Nb is the number of simulated B particles. If the collision is treated by the 
hard-sphere model, the scattering is isotropic. The post-collision velocities are then: 
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where the two R ’ s are the same unit vector with a random direction. However, the 
hard sphere model is generally insufficient for most applications because it does not 
match the gas’s viscosity or diffusion rate as a function of temperature. On the other 
hand, the most accurate Lennard-Jones potential is not suited for the DSMC method 



because the collision calculation is very time-consuming. But, the simple inverse- 
power potential V (r) = a / r a is sufficiently accurate and very efficient for collision 

calculations. The differential cross section o for the inverse-power potential takes the 
fonn: 


a = Ag- A ' a (2.2.7) 

where g is the relative speed and A is a function of the deflection angle X- If A is 
independent of x, the probability odil/o T that the post-collision relative velocity is 

directed into the solid angle dQ becomes <722/471, or the scattering is isotropic and the 
post-collision velocities are the same as that of the hard-sphere model. The variable 
hard-sphere (VHS) model was proposed by Bird [41] to avoid the divergence of the 
total collision cross section and using isotropic scattering for the inverse-power 
potential. In the VHS model, the coefficient A is constant and, once a gas is selected, 
A and a can be chosen to best fit the Chapman-Enskog viscosity to experimental data. 
This leads to a power law temperature dependence of the coefficient of viscosity such 
that: 
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(2.2.9) 


and co is bounded by 0.5 and 1.0 at the hard-sphere and Maxwell molecule limits, 
respectively. The VHS parameters used in the current study are discussed in Section 
5.3. 

2. 2. 1.1 Simple Gas 

For a single species gas in a computational cell, let ,v 2 . . ,v N be the velocities of 

the simulated particles. The probability P, that molecule i collides in At with some 
other particle is given by Equation (2.2.5). For a simple gas, this is: 

N N 

= S,Ms„)^ = Ip, (2.2.10) 
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where g tj = v ; . - Vj and 


P lJ = N-'ng IJ o r (g IJ )At 


(2.2.11) 


Note that P n = 0 because g„ = 0. The term for j = i should therefore be omitted in 
Equation (2.2.10). 



If gij is then replaced by g max , the resulting probability is: 


(n)„ = N " n 8™ a Ag m )*t (22.12) 

This probability is greater than Py. For Py = (Py) max , Equation (2.2. 10) takes the 
fonn: 


(^)_ =nN-'(N-l) g ^o T ( gmx )A, (2.2.13) 

The number of collisions in At for ( Pi) max is 

= (2.2.14) 

i = 1 

The factor of 14 on the right hand side of Equation (2.2.14) is to avoid double 
counting collisions between like particles. N max is called the maximum collision 
number. The computing time necessary to obtain the real collision number 

N L ,„ = itt p u < 2 - 2 - 15 > 

1=1 j=l 

is proportional to N 2 . This is not insignificant because collisions in all cells in a 
computational domain must be calculated at each time step. In the null-collision 
method, the task is avoided by regarding N max as the collision number in At. To 
compensate for excess collisions, N max -N c , some of the N max collisions are regarded as 
null collisions. We now define: 


< 2 - 216 > 

where q fj = g i} (J r (g.. )/g taa O T (g max ) - This means that the pair (/,/') first makes a 

tentative collision with a probability of (Pij)max- The tentative collision results in a 
real collision with a probability of qy and a null collision with a probability of 1 -qy. 
The number of tentative collisions is (. Py) ma x times the number of pairs N(N- 1)/2; the 
product is equal to N max . A tentative collision pair can be chosen randomly because 
(P ij)m ax is common to all (/,/). So, to summarize, we first obtain N max . In general, 

N max is not an integer. Let p = iV max - \_N max J , where |_^max J is the integral part of 

N ma . x- We then set A max = |_/V max J + 1 with a probability of p and At max = [/V max J with 
a probability of 1 -p. The following two steps are then repeated N max times: 

1 . Choose a pair (/,/) of i ± j randomly and obtain qy. 



2. Call a random number U. If U < q,j, the pair collides. If they do collide, the 
pre-collision velocities are replaced by the post-collision velocities given in 
Equation (2.2.6) where the particles A and B are of the same kind. 

2 .2 . 1 .2 Gas Mixture 

For the sake of simplicity, consider a binary mixture of species A and B. The types of 
collisions are A- A, B-B, and A-B. The A-B collision is a collision between unlike 
particles. The other two types are treated as in the gas of a simple gas above. Let Na 
and Ns be the numbers of simulated particles A and B in the cell. The probability that 
molecule A f collides with some molecule B in At is 


P A =n B go AB At (2.2.17) 

where g is the relative speed and of is the total cross section for the A-B collision. 
The average is for all simulated molecules B in the cell, so the probability takes the 
fonn: 
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where 
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If we again introduce the idea of the maximum collision number where gAi.ru is 
replaced by gf x , we have the upper bound of Pa. in'- 
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The maximum number of A-B collisions in At for P AB X is 
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Note that the factor of Vi is not necessary for A-B collisions. Substitution of Eq. 
(2.2.20) into Eq. (2.2.21) gives us: 

Nl =n t N iS “^“ (sl)Al (22.22) 


The number density n B is expressed as: 
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where V c is the cell volume and W is the weight. As discussed previously, one 
simulated particle represents W real particles. Presently, assume that all the particles 
in a cell have a common weight. The collision procedure in the previous section is 
then repeated for all species combinations A- A, A-B, and B-B using the updated 
definitions of number of collisions and probabilities for like and unlike collision pairs 
as required. 

2.2.2 Internal Degrees-of-Freedom 

A typical application of the DSMC method involves the calculation of many billions 
of intermolecular collisions, and this leads to the requirement for computationally 
simple molecular models. This has led to the development of “phenomenological” 
models that focus on the generation of accurate results for the observed gas properties 
but do not attempt to provide the most realistic model of individual collisions at the 
microscopic level[13]. The developers of classical kinetic theory faced a similar 
problem in that they required mathematically tractable molecular models that 
permitted the derivation of analytical expressions for the macroscopic transport 
properties. The classical models were the obvious choice for early DSMC work, but 
were found to have serious shortcomings that have been overcome by 
phenomenological models that have been introduced in the context of the DSMC 
method. 

We have already seen a phenomenological model be applied to elastic collisions by 
way of the variable hard sphere model, which combines the unrealistic but simple 
isotropic scattering of a hard sphere with a realistic variation of diameter with the 
relative speed of the collision pair. This model leads to accurate results for viscosity 
and heat conduction coefficients. 

In the same way, classical kinetic theory did not lead to any useful models for 
inelastic collisions. For example, the rough sphere model retains the deficiencies of 
the hard sphere model with regard to the transport properties, and it is unable to deal 
with the quantum effects that cause most gases to have fewer than three rotational 
degrees of freedom and has a fixed and unrealistically fast rotational relaxation time. 
This problem was solved in the context of the DSMC method by the introduction of 
the Larsen-Borgnakke, or L-B, model[42] for rotation which can be regarded as the 
archetypal phenomenological procedure. It is an add-on to the VHS model and, for a 
fraction of the collisions that are chosen to match measured rotational relaxation 
rates, the post-collision rotational energies are selected from the equilibrium 
distribution that correspond to the collision energy. The L-B model was extended to 
the vibrational modes, but it assumed that the internal energies were continuously 
distributed and was found lacking. The situation was corrected by the 
introduction[43] of the quantum version of the L-B model. 

DSMC methods test colliding particles for possible internal energy exchanges or 
chemical reactions. The collision selection rules use probabilities that must be 
defined such that they result in the desired macroscopic relaxation rate behavior. This 
typically involves identifying a macroscopic rate equation that the simulation 



attempts to match through appropriate definitions of selection probabilities P. These 
rate equations often contain parameters such as collision numbers Z, characteristic 
times x, or rate coefficients that are estimated from experiments or analysis. 
Ambiguity regarding specific definitions of the targeted macroscopic rates and 
coefficients has often led to confusion in defining the appropriate selection 
probabilities in DMSC methods. 

In studies of vibrational relaxation, Haas [44] developed a relationship between the 
appropriate DSMC selection probability P and the resulting Landau-Teller 
macroscopic relaxation behavior. This relationship was later generalized, allowing 
for rotational relaxation, by Lumpkin et al. [45] in the fonn of a “correction factor” for 
collision numbers Z. It was observed that the DSMC probabilities P, in general, 
cannot be specified by the assumption P=l/Z, which has been employed throughout 
the DSMC community. Dissimilar applications of the selection probability and 
incompatible definitions of collision number led to disagreement among DSMC 
researchers regarding these correction factors. 

Haas et al. [46] went on to clarify the definitions of macroscopic rate parameters 
employed in rotational relaxation, including differing interpretations of Jeans’ 
equation[47]. It was shown that relaxing energies, when plotted against the 
cumulative number of molecular collisions, lead to a single solution, regardless of the 
assumed intennolecular potential in an adiabatic bath. Differing definitions of 
collision time T c between experiments and DSMC were clarified. The two most 
widely used selection methodologies were discussed, and DSMC rotational relaxation 
probabilities P, specific to each methodology, were derived to match the macroscopic 
rate equations exactly. Prevalent selection methods cannot reproduce Jeans’ equation 
exactly for the case of heteromolecular collisions, nor can they reproduce Jeans’ and 
Landau-Teller[48] rate behavior for simultaneous rotational and vibrational 
relaxation. Thus, a new selection methodology was introduced that was better suited 
to the general case of gas mixtures experiencing multimode internal energy 
relaxation. Although constant Z values were used, the developments were valid for 
variable Z as well. 

Gimelshein et al.[ 49, 50] proposed a simplified inelastic collision algorithm using 
particle selection prohibiting multiple relaxations. They felt that the selection 
methodology proposed by Haas [46] lead to an unnecessarily complicated functional 
dependence between the probability and the collision number. To avoid the 
complexity, a straightforward selection methodology was developed which also 
resulted in relaxation behavior matching the Landau-Teller equation when the correct 
correction factor was used. 

Gimelshein et a/.[50] took a closer look at vibrational relaxation and derived an exact 
relationship that connects the vibrational relaxation number, Z vi b D M , used in the 
DSMC method and that employed in continuum simulations. An approximate 
expression for Z vi b DSMC was also derived that is cost-effective and applicable when 
translational temperature is larger than the vibrational temperature. 



2.2.2. 1 Rotational and Vibrational Energy Transfer 


Models for rotational and vibrational energy transfer with translational energy in a 
collision have been explored previously. The Larsen-Borgnakke (L-B) model[42] is 
usually used to describe the rotational and vibrational energy transfer in the DSMC 
method. A fraction of the total number of collisions, 1/Z VV /, or 1 IZ rot , lead to the 
energy exchange, where Z v ;& and Z rot are the vibrational and rotational collision 
numbers, respectively. Originally, the internal energy modes were treated as 
continuous. These energy modes are, of course, discrete. To address this issue, 
Bergemann and Boyd introduced the algorithms of the L-B model with discrete 
vibrational and rotational energy levels for diatomic molecules [43, 51]. Gimelshein 
et al. then extended the L-B model for the case of polyatomic gas mixtures with 
discrete vibrational and rotational energies[49]. 


The temperature dependent nature of rotational and vibrational relaxation rates has 
been developed by Parker[52] and Millikan and White[53], which provides the 
characteristic relaxation times as functions of the temperature for a wide range of 
pure and mixed diatomic gases. In the L-B model, the vibrational and rotational post- 
collision energies are modeled according to the local equilibrium distribution 
functions. The derivation of these models for application in the DSMC method can 
be viewed in Reference [2], but the vibrational relaxation model resulting from 
Millikan and White, which is used in this study, is: 
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where C/ and C 2 are constants used to fit the experimental data. 

However, these lead to the incorrect macroscopic relaxation rate. Lumpkin and Haas 
et al. [45, 46] introduced a correction factor for the rotational relaxation rate in the 
DSMC method. Most continuum analyses and computational codes use Jeans’ 
equation[47] to model rotational relaxation. In a single species gas, Jeans’ equation 
uses the instantaneous equilibrium rotational energy E rot (t) and a time constant, 
which can be written as the product of the instantaneous collision time x c and a 
characteristic collision number Z as: 
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where E rot is the mean rotational energy per molecule in the ensemble, and t 
represents time. Denoting the rotational degrees of freedom by ^ rot , E ro , is defined 
by: 
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and will change with time as the translational temperature T varies during relaxation 
under adiabatic conditions. Similarly, z c may vary with T, depending upon the 
intermolecular potential assumed for the gas. Most DSMC methods employ the VHS 
model, which assumes rigid-sphere collision mechanics, but establishes the collision 
rate in a manner consistent with an inverse-power potential. The relative translational 
energy of colliding particles near equilibrium, as a result of biasing due to collision 
selection, is distributed with ^,-ans degrees of freedom given by: 

4 --=5 -2® (2.2.27) 
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Note that this can differ from = 3 associated with the translational energy of an 
equilibrium ensemble of particles in a reservoir. Biasing from collision selection 
plays a significant role in the relaxation of the energy in a reservoir and becomes 
apparent in the derivation of the relaxation probability. Haas[46] derived the 
relationship for the particle selection method for a single species gas (allowing for 
multiple relaxation) as: 
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that correctly reproduced the relaxation rate given by Jeans’ equation. He goes on to 
say that P pan = 1/Z could be a fair assumption for adiabatic relaxation, but only if one 
redefined Z by using adiabatic relaxation data that has been cast into the isothermal 
fonn of Jeans’ equation as suggested by Bird [2], As such for adiabatic flows, P par t = 
(1 IZsird) = (l/Z)(l+£ rot /£r) can be a fair assumption to the exact expression. Better 
yet, the probabilities were shown to approach asymptotic limits described by: 
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which represents an excellent approximation when Z > 10 for the particle-selection 
methodology. 

Gimelshein[50] then analyzed the vibrational relaxation rates in the DSMC method. 
The main difference between rotation-translation and vibration-translation relaxation 
is that while the rotational mode can be describes as having a constant number of 
degrees-of-freedom at all temperatures except cryogenic, the number of vibrational 
degrees-of-freedom varies with the temperature. For a simple harmonic oscillator, 
the number of vibrational degrees-of-freedom ^ v ,b is: 
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where 0 v ib is the characteristic vibrational temperature of the mode. This means that 
the relationship for rotational relaxation probability given above are not applicable to 
the vibration-translation energy transfer. Gimelshein derives the relationships that 
allow one to match the vibrational relaxation rate in DSMC with that predicted by the 
Landau-Teller equation[48]: 
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where t is again time, E vi b is the mean vibrational energy per molecule in the 
ensemble, and E vi b is the instantaneous equilibrium vibrational energy, defined as: 


C (T ) 

17* vib \ trans / 

^ vib ~ ^ tn 


(2.2.32) 


He goes on to derive the expression for the DSMC collision number for vibration as: 
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The probability of a vibrationally inelastic collision in the DSMC method is then 
P = 1 / ZJ' . Although exact, Equation (2.2.33) is costly for direct use in DSMC 
computations because it involves the calculation of the vibrational and translational 
temperatures, and then the subsequent numerical detennination of T’ (post collisional 
temperature). However, the exact expression can be significantly simplified when the 
vibrational and translational temperatures are close: 
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The above equation is applicable only when the difference between the vibrational 
and translational temperatures is small. However, it was found that if the 
translational temperature is used as T, Equation (2.2.34) reproduced well the correct 
expression in Equation (2.2.33) for all cases where the vibrational temperature was 
lower than the translational one. 

2.2. 2.2 Electronic Energy Levels 

In the past, Bird[ 1 5] and Carlson et a/. [17, 25] modeled electronically excited levels 
and electronic transitions in DSMC. In their work, each excited particle was assigned 
a distribution over all available states, assuming that the electronic energy mode was 
in equilibrium. Gallis and Harvey[16, 54] modeled non-equilibrium thermal radiation 
in DSMC as well as electronic excitation processes. However, for their method, 
accurate cross sections for electronic excitation and absorption were required and the 
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uncertainty in the cross sections is not well known for all the species. Because of 
these deficiencies, a much more accessible non-equilibrium electronic energy model 
is developed in this thesis in Section 3.1. 

2.2. 2. 3 Selection Methodology 


Haas [46] derives the relationship for a versatile L-B selection methodology that is 
applicable to gas mixtures and is adapted from the particle-selection method; 
probability, P, is applied to each molecule individually within the colliding pair, but 
the event in which both particles relax is prohibited (called particle-selection 
prohibiting double relaxation). This ensures that relaxation of A molecules is 
independent of the energy of B molecules, and the L-B mechanics can therefore 
reproduce the relaxation rate behavior for each species independently. The 
appropriate functional form for the relaxation probability was defined as: 
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where the Z values are from continuum methods. This expression involves only 
rotational energy exchange. The expression for multimode energy relaxation is 
presented by Haas [46], and the reader is referred to the reference for the complete 
derivation. 


Gimelshein[49, 50] believed that, although correct, the above selection methodology 
leads to an unnecessarily complicated functional dependence between the probability 
and the collision number. This is especially true when multi-mode relaxation is 
considered. To avoid the complexity, a straightforward selection methodology was 
used that is illustrated in Figure 2.1. The first step is to select a random number Rn 
uniformly between 0 and 1 . Then, using the value of Rn, the possible relaxation 
events are tested, so that either rotational or vibrational relaxation may occur for 
particle A or particle B. Multiple relaxation events are prohibited. 


o 1 




^4 — ^\,A + Pf,B + ^ > w.A + ^ > v,B 

Figure 2. 1 : Selection methodology for inelastic and elastic collisions. 

Similar to that suggested by Haas, the simplified selection methodology allows only 
one relaxation event to take place in a collision. Vibrational and rotational relaxation 
rates of different species in the gas mixture are independent, and the necessary 
relaxation rates may also be accurately captured. It can be easily implemented and is 
characterized by a simpler dependence between the collision numbers and the 
relaxation probabilities. This method should not be applied when the total probability 
of internal-translational energy exchange is larger than unity. If the sum of the 
probabilities is greater than unity, an alternative methodology is implemented as 
shown in Figure 2.2. These selection methodologies are easily updated to include 
electronic energy exchange. 
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Figure 2.2: Selection methodology for inelastic collisions if the sum of the 
probabilities is greater than unity. 

2.2.3 Chemical Reactions 

Once the collision partners have been selected as described in Section 2.2.1, they are 
then tested to see if a chemical reaction will occur. A typical bimolecular reaction 
may be written as: 


A + B <->C + D 

As long as the reaction takes place in a single step with no species other than the 
reactants present, it has been found that the rate equation for the change of 
concentration of species A may be written: 

-^T = k f {T)n A n B -k r {T)n c n D (2.2.36) 

at 

The rate coefficients are functions of temperature and are independent of the number 
densities and time. A further result that originally had an empirical basis is that the 
rate coefficients are of the form: 
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where a and b are constants and E a is called the activation energy of the reaction. 

The Total Collision Energy (TCE) model[55] of Bird was developed in order to 
obtain a rate coefficient based on a chemical reaction cross section which takes the 
fonn: 
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where Ci and C 2 are constants, E co n is the collision energy available to the reaction, 
and £ is the average number of internal degrees of freedom that contribute to the 
collision energy. The resulting rate coefficient becomes: 
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where the constants Ci and C? are: 
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The probability of the reaction then occurring at each collision with E c > E a between 
power-law VHS particles is then: 
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The same can be done for tennolecular, dissociation/recombination reactions 
represented as: 


AB + T A + B + T 

where AB is the dissociating molecule and T is the third body particle. The rate of 
fonnation of either species A or B in the forward reaction can be written: 
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while the loss of species A in the reverse reaction is 

dn. 


dt 


A = k r n A n B n T 


(2.2.43) 


The corresponding reaction probability for the forward dissociation reaction is the 
same as that in Equation (2.2.41) and the probability of a recombination can be 
written: 
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where the activation energy and, therefore, contributing internal degrees of freedom 
are taken as zero. 

These are the typically implemented set of probabilities for both forward and reverse 
reactions in DSMC codes. They require available data for each forward and reverse 
reaction. However, these data are not always readily available for the reverse 
reactions. Also, many CFD codes make use of statistical mechanics in order to 
specify the reverse reaction rate. If comparisons are to be made to CFD results, it 
follows that the DSMC code should use statistical mechanics to obtain the reverse 
reaction rate where data are not available. This methodology is not used in general 
where data are available because the reaction probabilities are functions of the 
macroscopic temperature, which is less physically realistic than the dependence on 
the collision energy for the forward reactions. Another use of statistical mechanics is 
to calculate the reverse reaction rate from the forward reaction data, and then choose 
the relevant coefficients from the available reverse reaction data to most closely 
match the rate computed from statistical mechanics. The following modifications 
were implemented in the DAC software for those reverse reactions where data were 
unavailable. 

Chemical equilibrium requires a balance between the forward and reverse reactions 
such that the concentrations of all species remain constant with time. Therefore, for a 
bimolecular reaction, 
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where K eq is the equilibrium constant. Similarly, the equilibrium constant for the 
termolecular dissociation-recombination reaction is written as: 





Equilibrium statistical mechanics relates the numbers of the various species in a 
system to the partition functions. For the bimolecular reaction. 


N C N D 

n a n b 


C /~iD 


Q L Q 
Q A Q 1 


-exp 


E \ 

a_ 

kT , 


(2.2.47) 


where Q is the partition function for the species that is indicated by the superscript, 
and E a should be interpreted as the difference between the forward and reverse 
activation energies. The number N is related to the volume V of the system by N = 
nV, which results in the bimolecular equilibrium constant written as: 
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For the tennolecular reaction, the system volume does not cancel and the 
corresponding equation is: 
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The partition function for each species can be written as the product of the separate 
translational, rotational, vibrational, and electronic partition functions: 


Q = QtransQrotQvibQel (2-2-50) 

First, the translational partition function is: 

Q m = v(2xmkT/h 2 f (2.2.51) 

The rotational partition function is: 
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where e’ is a symmetry factor equal to 2 for homonuclear molecules and 1 for 
heteronuclear molecules. The harmonic oscillator model leads to a vibrational 
partition function of: 
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Finally, the electronic partition function is summed over the j levels as: 
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where g, are the degeneracies. Finally, it can be shown[2] that, for a bimolecular 
reverse reaction given the data for the forward reaction ( a ’s and />’ s), the reverse 
reaction probability is: 
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For the case where the reverse reaction has finite activation energy, the probability is 
then: 
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2.2.4 Electric Field 


Various researchers have modeled the electric field structure in a hypersonic shock 
layer in the past. The first such calculation was by Bird [26]. In his work, the 
electrons were moved through the computational grid with the ions that they were 
created with. No explicit electric field calculation was made, so the ions and 
electrons were not accelerated, and the simulation time step was that of the heavy 
particles. Later, Bird modified this method[27] to include a calculation of the 
ambipolar electric field using a form of the equation derived by Langmuir and 
Tonks[56] 
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In this equation, T e and n e are the average electron translational temperature and 
number density, macroscopic quantities derived from the DSMC solution. This 
equation can be derived from the macroscopic conservation of momentum for the 
electron species using the assumption of negligible inertial effects, negligible friction 
due to collisions, zero magnetic field and constant translational temperature. In his 
work, Bird calculated E am bipoiar using the results from a previously converged 
calculation. He then re-converged the DSMC calculation, applying the prescribed 
acceleration to the charged particles at each time step from E am bipoiar- The field 
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Eambipoiar was calculated from the re-converged DMSC solution, and the entire process 
repeated until there was no change in the resulting flow field parameters. The 
movement of the electrons was still tied to that of the ions in this method. 


In the work by Gallis and Harvey[28], the electrons and ions were again moved 
together, and the electric field was calculated in a similar manner to equation (2.2.58). 
Here, however, the assumption of isothennal electrons was not made: 
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and the gradient of the electron pressure was instead used to calculate the ambipolar 
electric field. Gallis and Harvey used a simulation time step corresponding to the cell 
crossing time of the heavy particles, and accelerated the heavy particles according to 
the field E am bipoiar ■ They then computed the average ion velocity in a given cell, and 
adjusted the average electron velocity in the same cell to match. In this way, the 
electron energies were affected by the electric field in an average sense. 

Carlson, Hassan and Taylor[25, 29] devised a method in which assumptions of zero 
net current and charge neutrality were used to calculate the electric field without the 
use of any macroscopic quantities. They wrote down Newton’s law for a charged 
particle of species s moving under the influence of an electric field: 

m s — = q E (2.2.60) 

dt 

and thus the expressions for average electron and ion velocity by summing over 
individual particles: 
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where E is the average electric field over a simulation time step, N is the number of 
simulator particles, v 0 is a velocity vector at the start of a simulation iteration, and the 

subscripts e and i refer to electron and ion species, respectively. By enforcing charge 
neutrality, N e = A), they then solved for the ambipolar electric field required to make 
the average ion and electron velocities equal, while enforcing the assumption of zero 
net current: 
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The electric field given above was used to accelerate the charged particles during the 
move portion of the DSMC algorithm. Additionally, a simple model of the particle 
acceleration in the plasma sheath at the vehicle surface was implemented. The 
electrons were reflected specularly from the vehicle surface to simulate their 
reflection in the strong negative potential gradient that would exist in the plasma 
sheath. An additional energy increment, e<p, was added to the surface heat transfer 
measurement for each ion impacting the surface. The expected potential drop across 
the sheath, (f), was calculated from one-dimensional collisionless sheath theory. 

Boyd [30] developed a model for the electric field of a weakly ionized plasma that is 
similar to the one originally proposed by Bird. The model was first implemented in a 
simulation of the plasma flow through an arcjet type thruster. Each electron is moved 
with the velocity of an ion throughout the domain, and the charged particles do not 
receive any velocity increment because of their response to an electric field. 

However, in this model, the average ion velocity of the cell in which the electron is 
located is used to move the electron particles, and it is not associated with a specific 
ion particle. This model maintains charge neutrality in an approximate sense, and 
results in lower computational overhead than the first model proposed by Bird. 

In order to calculate the self-induced electric field without any assumption of charge 
neutrality, zero net current, or ambipolar diffusion, the electrostatic Poisson equation 
must be solved. This was carried out for very rarefied flows surrounding spacecraft 
in low Earth orbit by two groups of researchers[57, 58]. While interesting, the 
preliminary results presented by these researchers in the cited references are for flow 
fields of significantly different structure than the shock layers being considered in this 
work. Gallis et a/.[ 59] presented preliminary results in which they solved the 
electrostatic Poisson equation in conjunction with the DSMC method for a hypersonic 
helium flow field. Unfortunately, the results presented in this work were sparse, and 
the physical processes associated with reacting air were not included in the analysis, 
and a systematic study of the effect of rigorously including the electric field in the 
computation on the flow field parameters was not presented. 

Presently, a hybrid of Bird’s original model and Boyd’s model[30] was implemented. 
The electrons were kept separate but retained their own velocities as was done in 
Boyd’s model; however, they are moved along with their associated ion as was done 
in Bird’s model. By doing this, all particles are moved at the heavy particle time 
scale, but collisions occur at the electron time scale. 

2.3 Recent Direct Simulation Monte Carlo Advances 

Molecular-level chemistry models that predict equilibrium and non-equilibrium 
reaction rates using only kinetic theory and fundamental molecular properties (i.e., no 
macroscopic reaction-rate information) have recently been proposed[ 10, 11, 13, 14]. 
In addition, a new vibrational relaxation model has also been introduced that only 
relies on fundamental molecular properties and a single measure of vibrational 
collision number at a reference temperature[10]. Preliminary evaluations [14] indicate 
that the resultant reaction rates are in very good agreement with measured Arrhenius 
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rates for near-equilibrium conditions and with both measured rates and other 
theoretical models for far-from-equilibrium conditions. 

2.3.1 Vibrational Relaxation 

A new vibrational relaxation model has recently been introduced by Bird[ 1 0] . The 
data for the vibrational collision number, Z vi b, is generally presented as a function of 
temperature in a Millikan-White plot. A convenient way of using this information in 
a DSMC procedure is through an expression that employs the values of Z vi b at two 
temperatures. The temperatures may be set to the characteristic temperature of 
vibration, 0 vi /„ and the characteristic temperature of dissociation, ©diss- Then, if the 
corresponding collision numbers are Z c and Z,/ m , the variation of the collision number 
with temperature T is given by: 
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The problem with the Millikan- White data for most gases is that it leads to collision 
numbers less than unity at the dissociation temperature. This is physically impossible 
and the most effective way to deal with the problem is to set the collision number to 
unity at the dissociation temperature. Also, it is preferable to set the second 
temperature, T re f, where the collision number is Z re f The collision number then 
becomes: 
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For a procedure based on this equation, the required data for each vibrational mode is 
the vibrational collision number at the reference temperature. The one-third-power 
law is preserved, but there is a decrease in the rate at higher temperatures. This is 
essentially a more precise alternative to the “Park correction factor ’[60] that is 
commonly employed in CFD. 

Consider a collision in which the sum of the relative translational energy and the pre- 
collision vibrational energy of the molecule under consideration is E co n. The 
“quantized collision temperature” is defined by: 
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where the collision energy has been quantized through the truncation of i max to an 
integer. The relevant parameters for the new vibrational model used are listed in 
Table A. 3. 

2.3.2 Chemical Reactions - The Quantum-Kinetic Chemistry Model 

Bird has proposed[10, 11, 13] a set of molecular-level chemistry models based solely 
on fundamental properties of the two colliding particles: total collision energy, 
quantized vibrational levels, and molecular dissociation energies. These models li nk 
chemical-reaction cross-sections to the energy-exchange process and the probability 
of transition between vibrational energy states. The Larsen-Borgnakke[42] 
procedures and the principle of microscopic reversibility are then used to derive 
simple models for recombination reactions and for the reverse (exothennic) reactions. 
These models do not require any macroscopic data, and they balance the fluxes into 
and out of each state, thus satisfying microscopic reversibility. As mentioned earlier, 
preliminary evaluations[14] indicate that the resultant reaction rates are in very good 
agreement with measured Arrhenius rates for near-equilibrium conditions and with 
both measured rates and other theoretical models for far- from equilibrium conditions. 

2.3.2. 1 Dissociation Reactions 

Bird suggests[10] that a complete and consistent description of vibrational energy 
exchange between colliding molecules must include dissociation. This concept is 
embodied in the energy-threshold model of Bird[2], According to this model, during 
the energy-exchange process of a diatomic molecule, if the energy content of the 
vibrational mode exceeds the energy threshold of the dissociation or exchange 
reaction, the reaction AB + T — > A + B + T occurs. 

In his critique of Bird’s model, Lord[61] pointed out that a “continuum” of states (i.e., 
an infinite number) exists above the dissociation limit. Therefore, Lord argues, if a 
dissociation reaction is energetically possible, it occurs. Thus, assuming a particular 
collision between two simulators, where at least one of them is a diatomic molecule, 
the serial application of the Larsen-Borgnakke energy exchange model would make 
the collision energy: 
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available to the vibrational mode of the molecule in question. In this expression of 
collision energy, E tram , pa ir is the translational energy of the pair, and E vib mol is the 
vibrational energy of the molecule considered for dissociation. Assuming a hannonic 
oscillator model, the maximum vibrational level that can be reached is that given in 
Equation (2.3.3), where the vibrational collision number has been defined as unity. If 
this level is higher than the dissociation level (Q = 0/ / 0,,), or 
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a dissociation occurs. 


When the serial application of the quantum Larsen-Borgnakke model is considered, 
potential post-collision states i are chosen uniformly from the states equal to or 
below i max and are selected through an acceptance-rejection routine[2] with 
probability 
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where co is the temperature exponent of the coefficient of viscosity. Then, the 
corresponding rate coefficient in an equilibrium VHS gas at temperature T is: 

kf (T) = gR“fr(i aa ) A ‘ T (2.3.7) 


The degeneracy g of the reaction is included to cover the case of a polyatomic 
molecule in which the reaction can occur through more than one mode with the same 
characteristic vibrational temperature. The collision rate parameter R^ u ' is the 
collision rate for collisions between gas species AB and T divided by the number 
density product. It can be written for an equilibrium VHS gas as 
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where r re f is the molecular radius at temperature T re f and m r is the reduced mass. For 
dissociations, it is necessary to add the symmetry factor e that is 1 for homonuclear 
and 2 for heteronuclear molecules. The parameter T (/ max ) ABJ is the fraction of 

collisions between AB and T that have sufficient energy to meet the condition of 
equation (2.3.5). With i as the pre-collision vibrational state of AB, the result for an 
equilibrium VHS gas is 
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Here, Q(a,x) = r(a,x)/T(a) is a fonn of the incomplete Gamma function and 
7,n)=i/[i - exp(-0 v /T )] is the vibrational partition function. As an example, the 

dissociation rate of molecular nitrogen is presented in Figure 2.3 for the Q-K model 
as compared to several rates from the literature, the first of which is the value used in 
the DSMC simulations using the TCE chemistry model. 
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Figure 2.3: Comparison of Q-K dissociation rate of the reaction N 2 + N 2 — > N + N + 
N 2 with those from literature. 

2. 3. 2.2 Recombination Reactions 

The condition for recombination in a collision between species A and B is that 
another atom or molecule is within the “ternary collision volume” V co u. The 
probability of recombination is therefore 

Prec = nV coll (2.3.10) 

where n is the number density. The corresponding rate equation is 

K(T) = R>V CC „ (2.3.11) 

The ratio k f (T)/ k r ( T ) from Equations (2.3.7) and (2.3.1 1) may be equated to the 

equilibrium constant from statistical mechanics in order to determine the ternary 
collision volume. It has been found that a volume in the form 

V coll = aT h V ref (2.3.12) 

allows nearly an exact fit of this ratio to the equilibrium constant. The coefficients 
used for these adjustments can be viewed in Table A. 16. The reference volume V re f is 
most conveniently set to the volume of the sphere with radius equal to the sum of the 
reference radii of the three particles. 

The values of the constants a and b are determined in the context of the macroscopic 
temperature, but the evaluation of the ternary collision volume in the DSMC 
implementation is based on the quantities associated with the collision. The 
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“collision temperature” based on the relative translational energy in a collision with 
relative speed g is: 
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As an example, the cumulative recombination rate of two atomic nitrogen atoms is 
presented in Figure 2.4 as compared to values from the literature, the first of which is 
the rate computed by statistical mechanics from the dissociation of molecular 
nitrogen. 
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Figure 2.4: Comparison of Q-K recombination rate of atomic nitrogen with those 
from literature. 

2. 3. 2. 3 Exchange Reactions 

Exchange reactions are binary reactions written as A + B^C + D in which A and C 
are the molecules that split, while B and D are the atoms. The phenomenological 
DSMC procedure for this reaction is analogous to that for the dissociation model. For 
both the forward and reverse reactions, the reaction probability is equal to that of 
Equation (2.3.6) with i set to the vibrational level corresponding to the activation 
energy of the reaction, i a 
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This is a relative probability that must be normalized through division by the sum of 
probabilities from the ground state to i max , resulting in 
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This is sufficient if the level i a is much larger than unity, but some endothennic 
reactions and all exothermic reactions have either zero or small activation energies. 
The use of discrete levels then leads to unacceptably abrupt changes in the rate 
coefficients and it is impossible to obtain a near exact agreement between the ratio of 
forward to reverse rate coefficients and the equilibrium constant. The level i a is the 
highest level with energy below the activation energy E a and a smooth curve is 
obtained for the rate coefficients if the energy i a kO v in the numerator is replaced by 
E a . The DSMC probability of an exchange reaction in a collision with E c > E a is then 
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An exact expression can be written down for the fraction of collisions that satisfy the 
condition of Equation (2.3. 16) in an equilibrium VHS gas. With j as the pre-collision 
vibrational level, this is 
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and the corresponding rate coefficient is 

KT) = gR„u(dNIN) (2.3.18) 

This is a complex expression not easily programmed, but it is possible to develop an 
alternative expression that is based on physical reasoning rather than fonnal kinetic 
theory. The selection of possible post-collision vibrational levels in Equation (2.3.6) 
is from an equilibrium distribution so that the rate coefficient is simply the product of 
the degeneracy g, the collision rate parameter, and the fraction of the vibrational level 
i a in the Boltzmann distribution. The term i a k0 v in the exponential is replaced by E a 
in order to maintain consistency with Equation (2.3. 16), and the normalizing factor is 
simply the vibrational partition function, thus 
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The default activation energy for the forward reaction is the heat of reaction and that 
for the reverse reaction is zero. The ratio of the forward to the reverse rate is 
compared with the equilibrium constant and the two are brought into the best possible 
agreement by the upward adjustment of either the forward or reverse activation 
energy. These adjustments are temperature dependent through a relation similar to 
that in Equation (2.3. 12) for the ternary collision factor: 
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For the endothennic reaction, E exc h is the negative of the activation energy, which 
would then result in the second part of the equation adding the prescribed energy to 
the activation energy. For the exothennic reaction, E exc h is equal to the activation 
energy of the endothennic reaction, so the second part of the equation above does not 
come into effect. The coefficients used in these adjustments are presented in Table 
A. 17. In addition to these phenomenological adjustments, there may be experimental 
or theoretical evidence for activation energies above the default values. Should a 
phenomenological adjustment fall short of the desired value, both the forward and 
reverse activation energies may be adjusted to keep the forward to reverse rate ratio in 
agreement with the equilibrium constant. As an example, the endothennic and 
exothermic reaction rates for the exchange reaction N 2 + O — > NO + N are compared 
to rates from the literature in Figure 2.5 and Figure 2.6, respectively. The Q-K results 
presented in these figures are in good agreement with those presented from the 
literature. 
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Figure 2.5: Comparison of Q-K exchange rate for the endothermic reaction Nt + O 
NO + N with those from literature. 
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Figure 2.6: Comparison of Q-K exchange rate for the exothermic reaction NO + N 
N 2 + O with those from literature. 
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Chapter 3 : Extending the Quantum-Kinetic Chemistry Model 

For missions where the shock-layer temperature reaches tens-of-thousands of degrees 
Kelvin, the importance of including electronic energy levels and ionization reactions 
becomes more pronounced. Although the treatment of electronic energy levels and 
ionization using the DSMC method is not a new development, most of these methods 
are based on measured, equilibrium rates, which are always questionable at high 
temperature, especially when applied to non-equilibrium problems. 

In Chapter 2, the Quantum-Kinetic Chemistry Model of Bird was introduced. One of 
the attractive properties of this model is the fact that it does not rely on any 
macroscopic information, leveraging instead off of kinetic theory and physical data 
corresponding to the molecules undergoing each individual collision. The purpose of 
this chapter is to introduce new methods of treating electronic energy level transitions 
and chemical reactions involving charged species following the Q-K methodologies 
utilizing not only the vibrational energy levels, but also the electronic energy levels of 
the atoms and molecules. 


3.1 Detailed Electronic Energy Level Model 

To model the distribution of electronic energy in the DSMC technique, three 
procedures must be defined. First, when a particle is introduced into a simulation, it 
is necessary to obtain a new electronic energy from a given distribution through 
statistical sampling. A new electronic energy is also required when a particle changes 
its electronic energy after a collision with a surface. Under equilibrium conditions, 
the distribution has the well-known Boltzmann form. Second, it is necessary to 
statistically sample a new electronic energy following a collision, which involves 
electronic energy transfer. Third, a method to reproduce the electronic energy 
transition rate must be defined. Detailed examples of applying this model at 
hypersonic (re)entry conditions in air are presented subsequently in Section 4.3. 

3.1.1 Equilibrium Sampling 


Each electronic energy level j has a distinct energy, Ej, and degeneracy, g,. The 
Boltzmann distribution for the electronic energy levels at a given temperature T gives 
the following result for the fraction of particles in level j: 


felU) 


N j= gjCxp(-Ej/kT) 

AT 'max 

£s,exp(-£,./kr) 


(3.1.1) 


This distribution is used when either creating a new particle or after a surface 
collision in the DSMC simulation at a boundary specified at the temperature T. 
However, it is not possible to sample an electronic energy level j directly from the 
distribution. Therefore, an acceptance-rejection procedure is employed by selecting 
values for j from the following distribution: 



( ' g I ^y(-EJkT) 

gjOxp(-Ej/k T ) 


(3.1.2) 


where J is the value of j for which Equation (3 . 1 . 1 ) is a maximum. Unfortunately, 
since the degeneracy of each level is a variable specific to each species, there is no 
way to detennine the maximum level a priori, so the maximum level must be 
searched for at each implementation or saved and recalled for a constant boundary 
temperature. The sampling of a new electronic energy level then proceeds as follows: 

1 . select at random an electronic energy level evenly distributed between 0 and 
Jmax where J max is the maximum possible energy level; 

2 . determine the value J for which Equation (3 . 1 . 1 ) is a maximum; 

3. accept the value of j if f e i ’(j) > RAND, a random number; 

4. if the value of j is not accepted, return to step 1 . 

3.1.2 Post-Collision Sampling 

A phenomenological approach is usually adopted in the DSMC method when a 
collision occurs that involves energy transfer. The Larsen-Borgnakke method[42] 
samples a post-collision state from a combined distribution of the translational and 
electronic collision energies of the colliding particles. Based on the approach of 
Bergemann and Boyd[43] and Boyd[51], the Dirac delta function is employed to 
define the distribution of energies in Equation (3.1.1) in the following continuous 
form: 
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Consideration must now be given to the distribution of translational energy of the 
colliding particles. This distribution is naturally affected by the intermolecular model 
used since this determines the collision probability. For the present study, the 
variable hard sphere (VHS) collision model of Bird[41] is used. However, it is 
straightforward to develop the formulation for an alternative collision model such as 
the variable soft sphere of Koura[62]. For the VHS model, the distribution of the 
translational energies is: 
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Using Equations (3.1.3) and (3.1.4), the combined distribution for sampling a post- 
collision electronic energy level j ’ from the total collision energy E co n = E trans +E e i = 
Etrans ’+E e /’ is: 
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In applying the general Larsen-Borgnakke scheme, it is assumed that local 
thermodynamic equilibrium prevails. Therefore, the temperature T in Equation 
(3.1.5) is constant. Also, the total collision energy is constant, so it is only necessary 
to perform sampling of the post-collision state from the following distribution form: 

g(f\E^g r (E M -Eif (3.1.6) 


Again, an acceptance-rejection procedure is used. The normalized distribution that is 
required is obtained by finding the value of j’ for which Equation (3.1.6) is 
maximum, J\ which is different for each value of E co n. Therefore, the following 
distribution is obtained: 
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The procedures are implemented as follows: 

1 . given a pair of particles with total collision energy E co u which undergoes 
electronic energy exchange, determine/’; 

2. determine the maximum allowable electronic energy level obtainable from 

Ecoll, J”\ 

3 . take J* to be the smaller of J’ and J’ 


4. as described in the previous section on equilibrium sampling, perform an 
acceptance-rejection procedure to sample j ’ from: 
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3.1.3 Electronic Energy Level T ransitions 

A transition from electronic energy level i to level j for species A can be written in 
the form of a chemical reaction as A' + B — > A 1 + B . Following the Q-K model, the 
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simplest model is to assume that the transition occurs if the electronic energy level of 
A after a trial Larsen-Borgnakke redistribution of the collision energy, E co u = 
Etmns + E e i, is j = /. Since the transition is treated as a chemical reaction, there are 
J max transitions to consider where J max is the maximum energy level energetically 
possible. 

In the DSMC code, the process for determining if an energy level transition occurs, 
and if so, to which transition level, is as follows: 

1 . determine the value of J max and compute the denominator of Equation (3 .1.7); 

2. perfonn acceptance-rejection according to Equation (3.1 .7) until a test level is 
accepted. 

This selection process is analogous to the selection process used in selecting a 
chemical reaction from a list of possible reactions. 

3.2 Chemical Reactions Involving Charged Species 

As outlined in Section 2.3.2, Bird[10, 11, 13], supported by Gallis[ 14], proposes a set 
of molecular-level chemistry models based solely on fundamental properties of the 
two colliding molecules: total collision energy, quantized vibrational energy levels, 
and molecular dissociation energies. These models link chemical cross-sections to 
the energy exchange process and the probability of transition between vibrational 
energy states. The Larsen-Borgnakke procedures and the principle of microscopic 
reversibility are then used to derive simple models. The models do not require any 
macroscopic data and they function by seeking to balance fluxes into and out of each 
state, thus satisfying microscopic reversibility. Since reactions involving only neutral 
species have been covered thus far, processes involving charged species are presented 
in the current section. Again, more detailed comparisons for these reactions 
pertaining to hypersonic (re)entry in are be presented in Sections 4.4 through 4.7. 

3.2.1 Ionization Reactions 

The procedure for ionization is analogous to that for dissociation. The difference is 
that ionization is linked to the electronic energy level instead of the vibrational energy 
level as with dissociation. For ionization, the collision energy is defined as 
E co ii=E tra ns,pair+E e i e , pa rti C i e where the first term on the right is the relative translational 
energy of the colliding pair and the second tenn is the electronic energy of the 
particle being considered for ionization. If the collision energy is greater than the 
ionization energy of the particle, ionization occurs. The corresponding rate 
coefficient in an equilibrium VHS gas at temperature T is: 
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where the collision rate parameter is identical to that of Equation (2.3.8). The 
fraction of collisions between A and B that have sufficient energy to meet the 
requirement of ionization is then: 
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where z e i(T) A is the electronic partition function for particle A. 
3.2.2 Associative Ionization Reactions 


The forward (endothermic) associative ionization reaction can be thought of 
occurring as A + B — > AB — > AB ' + e where AB is an intennediate, excited 
molecule that may lead to the ionized species AB + plus an electron. The first step is 
similar to the recombination procedure described in Ref. [10]. If the AB molecule is 
in the ground vibrational state after a trial Larsen-Borgnakke redistribution of the 
collision energy, the intermediate molecule is “formed.” In the second step, the 
dissociation energy is first added to the original collision energy. The new collision 
energy is then distributed among the internal energy modes of the AB * molecule and 
then this molecule is tested for ionization as described above. If the remaining 
collision energy after redistribution plus the electronic energy of the molecule is 
greater than the ionization energy of the molecule, the associative ionization reaction 
occurs. The corresponding rate coefficient can then be written as: 

k?=l$S*%,/z»{Tf r (3.2.3) 


where the collision rate parameter is again identical to that of Equation (2.3.8). 
The vibrational partition function is included to account for the probability that the 
AB* molecule will be in the ground state. The fraction of collisions with sufficient 
energy for the molecule AB to ionize are then calculated by summing over all 
possible initial internal energy states of atoms A and B (electronic energy states) and 
all possible vibrational, rotational, and electronic energy states of the intermediate 
AB * molecule as follows: 
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The sums over the atomic electronic energy levels is over all possible levels. For the 
molecular internal energies, the maximum level is that which is the highest 
energetically possible. To start, an initial translational energy is assumed to be the 
average relative translational energy, or ( 5 / 2 — (o)kT . This is an approximation and 



may result in a discrepancy between the above analytical expression and the actual 
sampled rate. It is possible to integrate over all possible relative translational 
energies, but it makes an already long calculation (summing over all possible internal 
states) much longer. As will be shown subsequently, this formulation results in an 
acceptable comparison with the sampled values in Section 4.4. The first two 
expressions after the summations are simply the probabilities of atoms A and B being 
in electronic levels i e uf and i e i e B ', respectively. The next three terms are analogous to 
Equation (3.1.8) in that they are the probabilities of the respective internal mode of 
energy being at the specified quantum level at the specified collision energy. The 
collision energies are defined as: 
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where each successive collision energy subtracts the previous internal energy mode. 
Then the last term, as in the dissociation and ionization rate equations, is the fraction 
of collisions for which the translational energy is greater than the difference between 
the ionization energy of AB and the electronic energy of AB . 

The reverse (exothermic) associative ionization reaction is similar to the forward 
reaction. If the AB molecule is in the ground electronic state after a trial Larsen- 
Borgnakke redistribution of the collision energy, an intermediate molecule is 
“formed.” In the second step, the ionization energy is first added to the original 
collision energy. The new collision energy is again distributed among the internal 
energy modes of the AB molecule, which is then tested for dissociation as described 
above. If the remaining collision energy after redistribution plus the vibrational 
energy of the molecule is greater than the dissociation energy of the molecule, the 
reverse associative ionization reaction occurs. The corresponding rate coefficient can 
then be written as: 
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where the collision rate parameter is again identical to that of Equation (2.3.8). 
The last term is included to account for the probability that the AB * molecule will be 
in the ground electronic state. The fraction of collisions with sufficient energy for the 
molecule AB to ionize are then calculated by summing over all possible initial 
internal energy states of AB + (the electron has none) and all possible vibrational, 
rotational, and electronic energy states of the intermediate AB * molecule as follows: 
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The sums over the internal energy levels of AB' are over all possible levels. For the 
internal energies of AB , the maximum level is that which is the highest energetically 
possible. In order to do this, an initial translational energy is again assumed to be the 
average relative translational energy, or ( 5 / 2 — a))kT , along with the aforementioned 

cautions. The first three expressions after the summations are simply the probabilities 
of the AB 1 molecule being in the respective states. The next three terms are again 
analogous to Equation (3.1.8) in that they are the probabilities of the respective 
internal mode of energy being at the specified quantum level at the specified collision 
energy. The collision energies are defined as: 
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where each successive collision energy subtracts the previous internal energy mode. 
Then the last term, as in the dissociation and ionization rate equations, is the fraction 
of collisions for which the translational energy is greater than the difference between 
the dissociation energy of AB and the vibrational energy of AB 


3.2.3 Charge Exchange Reactions 


The charge exchange reaction is analogous to the exchange reaction outlined in 
Section 2. 3.2. 3. The binary exchange reaction can be written as A + B + A + + B . 
This reaction takes place when the electronic energy level of the colliding neutral 
particle, A, after a trial Larsen-Borgnakke redistribution of the sum of the relative 
translational energy and electronic energy of the neutral particle is equal to the 
electronic energy level with energy just below the activation energy of the reaction. 
The probability of the reaction occurring is then: 
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The corresponding rate equation for the forward reaction is then: 
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and for the reverse reaction is: 
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As with the exchange reactions in Section 2. 3.2. 3, the ratio of the forward and reverse 
reactions can be compared with the equilibrium constant and adjustments can be 
made to the phenomenological model to bring the model in closer agreement with 
statistical mechanics. 


3.2.4 Other Exchange Reactions 

In other exchange reactions involving charged species, the predominant mechanism is 
the breaking of a bond of one species (either neutral or charged) in order to exchange 
one atom to bond with the original atom, just as with an exchange reaction involving 
only neutral particles. Therefore, the logic for performing this reaction involving 
charged species is similar to that involving only neutral species. 

3.3 Selecting a Reaction 

Next, the Quantum-Kinetic chemistry model requires the selection of which reaction 
occurs. For many collision pairs, there may be more than one possible reaction. In 
the TCE chemistry model, the reaction probability is computed for each of the 
possible reactions and summed. This sum is then compared to a random number 
between zero and one. If the random number is less than the sum of the probabilities, 
a reaction is to take place. The selection of the reaction then proceeds by normalizing 
the reaction probabilities, and another random number is chosen to select which 
reaction is allowed. In the Q-K method, there is no clear analogy, and previous 
implementations have utilized a hierarchical approach where the reactions are 
considered by type first. For the neutral reactions, the dissociation reaction was 
always considered first and allowed to proceed even if an exchange reaction was a 
possibility. However, it was found that this approach tended to be biased in favor of 
dissociation reactions, especially at high temperatures where most of the collisions 
have larger relative translational energies and the vibrational level distribution is 
shifted to the higher levels. This method becomes even more undesirable when 
charged species reactions are included in the reaction set. Specifically, what happens 
when both dissociation and ionization reactions occur simultaneously? This event 
doesn’t actually occur for the current reaction set under consideration (only electron 
impact ionization of atomic species is considered herein), but this is obviously an 
issue that needs to be resolved. 
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An alternate approach to the hierarchical method was therefore implemented. This 
approach is, like the TCE model, based on assigning reaction probabilities to all of 
the possible reactions, summing them, then comparing the sum to a random number. 
If the random number is less than the sum of the probabilities, then the probabilities 
are nonnalized and then another random number is generated to select the reaction. 
The reaction probabilities for exchange and charge exchange reactions have already 
been given in Equations (2.3. 16) and (3.2. 10), respectively. For the remaining 
reactions that occur when the collision energy is such that the reaction occurs 
(dissociation, ionization, and associative ionization), the probability is set equal to 
one if it is determined that the reaction is allowed to occur. These probabilities are 
then summed and compared to a random number between zero and one (obviously, if 
a dissociation, ionization or associative ionization reaction are possible, a reaction 
will occur). If a reaction is to occur, P r i is defined as the reaction probability for the 
i th reaction, and the cumulative distribution for all reactions up to and including the j th 
reaction is: 





(33.1) 


where NR is the total number of reactions possible with the current reactants. A 
random number is then chosen to determine which reaction proceeds. 

As shown subsequently in Chapter 4, this selection methodology allows the sampled 
reaction rate to accurately reproduce the analytical rates shown in this chapter. There 
is one item that requires mentioning at this point that the chemistry set utilized in this 
thesis does not include. That is the possibility of a dissociation and ionization 
reaction occurring for the same set of particles. For example, when molecular 
nitrogen collides with an electron, the molecule can either dissociate or ionize. Under 
the rules for these processes listed above, when the respective criterion is met for 
either reaction, the probability of that reaction occurring is unity. What happens 
when both have a probability of unity? Using the above reaction selection 
methodology, it is shown in Figure 3.1 that the sampled rates still accurately 
reproduce the analytical rates. This occurs because the energy required for 
dissociation is less than that for ionization, so will always occur more frequently as 
shown in the figure. The highest reaction rate physically possible would be equal to 
the collision rate of the particles. As the limit where all of the collisions have 
sufficient energy for both reactions is approached, a departure from the analytical 
rates occurs since both are equal, resulting in sampling rates one half the analytic 
rates. 




Temperature (K) 

Figure 3.1: Sampled and analytic reaction rates when dissociation competes with 
ionization. 
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Chapter 4: Application of New Models to Upper- Atmosphere 

Hypersonic Flows 

In this chapter, the new electronic energy level transition and chemistry models are 
validated by examining their ability to reproduce rates of chemical reactions in the 
literature[5-7, 9, 15, 29, 60, 63-78]. The DSMC models are applied in a two- 
dimensional code where the results are accumulated over all cells to approximate a 
zero-dimensional simulation. The test gas is comprised of eleven species: O, N, CT, 
N 2 , NO, 0 + , N + , 0 2 + , N 2 + , NO+, and e'. The simulations were run at various 
temperatures between 10,000 K and 60,000 K and at a continuum number density of 
1.0e23 /nr’ in an adiabatic box 0.002 m on a side with 2 million particles. During the 
simulations, relaxation of the particles (including rotational, vibrational and 
electronic) was allowed to proceed as usual, but when a reaction was determined to 
take place, the number of reactions was advanced by one but the simulators were left 
unchanged, to maintain the thermodynamic energy balance with energy neither added 
nor removed from the gas. 

4,1 Measurement of Reaction/Transition Rates 

For a bi-molecular reaction A + B <-> C + D , the forward reaction rate is calculated 
as: 


dn, / 
dl 
n A n B 


(4.1.1) 


where the numerator is the change in number density of species A due to the reaction 
in the sampled time and the denominator is the product of the number densities of 
species A and B. The corresponding reverse reaction rate is: 


dn r / 

. _ / dt 

n c n D 


(4.1.2) 


For a tennolecular reaction AB + T <-> A + B + T , the forward reaction rate is 
calculated as: 




n AB n T 


(4.1.3) 


as the bi-molecular reaction rate is, but the reverse reaction rate is now: 



The capability of these techniques to reproduce measured equilibrium reaction rates 
will now be investigated. However, the uncertainty of the currently available 
measured rates for atmospheric hypersonic flows can be large. Most of the 
experimental data on which reaction rate sets are based were obtained for low 
temperatures (2,000 - 7,000 K). Even though uncertainty may be large, Park[60] 
suggests that the experimentally obtained reaction rates can be extrapolated up to 
30,000 K. 

The uncertainty of these measured rates is not always reported and in some cases is 
unknown. The problem becomes even more severe in the extrapolated portion of the 
temperature range (T> 10,000 K). Park[60] reports that some rates measured for the 
same reaction differ by more than one order of magnitude, which will be apparent in 
the figures that follow. Since this work is focused on the extension of the Q-K model 
to electronic energy levels and reactions involving charged particles, only these 
transition and reaction rates will be discussed in this chapter. However, comparisons 
of the analytical Q-K rate expressions, the sampled Q-K rates, and rates found 
throughout the literature for every reaction included in this study are presented in 
Appendix C. 

4.2 Definition of Translational, Rotational, Vibrational, and Overall 
Temperatures 

There are various ways to define temperatures with regards to the DSMC 
methodology. In this work, we will refer to statistical mechanics for the definition of 
temperatures by defining the energy as: 

£ = y 2 £RT (4.2.1) 

where e is the specific energy of the quantity of interest, £ is the number of degrees- 
of-freedom, R is the universal gas constant, and T is the temperature. By finding the 
energy, recognizing that k = mR , and rearranging Equation (4.2.1), one has the 
temperature of interest. For the translational temperature, we define £ tm ns as: 

(4-2.2) 

where c' 2 is the average square of the thermal velocity defined as: 

c 7I = ?-2 dc 0 + c 0 2 (4.2.3) 

where c is the mean molecular velocity and c 0 is the macroscopic stream or mean 
velocity. The translational temperature can then be defined as: 

r„ = ^ (4.2.4) 

since there are 3 translational degrees of freedom. 
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For the internal degrees of freedom of a specific species n, we define the average 
internal energy of a single particle as: 


F" - m f" 

^int rn nr / int 


and the temperature is then: 


2 K 

K’L 


(4.2.5) 


(4.2.6) 


For the overall internal temperature averaged over all of the N species, we then 
define: 


N 

y c" t 

/ 4 int li 


int N 

y r 

/ j bint 

n = 1 


(4.2.7) 


For rotational energy, the average energy of each species is found by summing the 
rotational energies of the particles in the computational cell and dividing by the 
number of particles. For the rotational degrees of freedom, it is assumed that the 
mode is fully excited and is set to 2 for diatomic molecules. 

For vibrational energy, the average energy of each species is again found by summing 
the vibrational energies of the particles and dividing by the number of particles. To 
then find the number of vibrational degrees of freedom for the species, using the 

definition of the average vibrational energy level i n = E" ib /kO vib and some 
manipulation, one gets: 

C, = 2i>(l + l/<“) (4.2.8) 


The electronic temperature is treated in a similar manner, but some modifications 
must be made to account for the discretized spacing of the electronic energy levels. 
To begin, the species electronic temperature is defined as: 
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(4.2.9) 


where the 0 and 1 denote the ground and first excited electronic energy levels, g is the 
degeneracy of the level, and n can be taken as the number density or total number of 
samples (result is the same since the ratio is being taken). This definition comes 
directly from the ratio of the Boltzmann distribution of the ground and first excited 



electronic energy levels. Then, the effective number of degrees of freedom is defined 
as: 


C n 

ele 


2 Ke 
kT 


ele 


(4.2.10) 


as we have defined the effective degrees of freedom of the other internal modes. 

4.3 Electronic Energy Level Transitions and Comparison to Equilibrium 

The following section applies the methodologies developed in Section 3.1 to the 
aforementioned equilibrium heat bath. There are excitation and quenching processes 
between each of the atomic/molecular energy levels, resulting in a large number of 
rates being required for previous electronic energy level transition models, both for 
CFD and DSMC. The new Q-K based transition model negates the necessity for each 
of these transition rates (most of which use approximate analytic formulas). In the 
following sub-sections, only a portion of the transitions presented in Appendix C are 
discussed for exemplary purposes only. 

4.3.1 Electron-Impact Excitation 

The process of electron-impact excitation is significant for both atoms and molecules 
in hypersonic shock-layer applications. To begin, the electron-impact transitions of 
atomic nitrogen are examined. The sampled transition rates are compared to values 
reported by Tayal[79], Park[60], and Frost[80]. The transitions from the first to the 
fifth level and from the second to the fifth level are presented in Figure 4.1 and Figure 
4.2, respectively. It is important to note that the 1-5 transition is considered a 
“forbidden” transition, while the 2-5 transition is an “allowed” transition. Even so, 
the phenomenological Q-K model compares very well to the reported rates from the 
literature. 
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Figure 4. 1 : Comparison of sampled and reported electronic energy level transition 
rates from the first to fifth energy levels of atomic nitrogen. 



Temperature (K) 

Figure 4.2: Comparison of sampled and reported electronic energy level transition 
rates from the second to the fifth energy levels of atomic nitrogen. 

For the electron-impact transitions of atomic oxygen, the sampled transition rates are 
compared to values reported by Zatsarinny[81], Park[60], and Bhatia[82]. Once 
again, even though one transition is allowed (level 1-5, Figure 4.3) and one is 
forbidden (level 2-5, Figure 4.4), the comparisons are quite favorable. 
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Figure 4.3: Comparison of sampled and reported electronic energy level transition 
rates from the first to fifth energy levels of atomic oxygen. 
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Figure 4.4: Comparison of sampled and reported electronic energy level transition 
rates from the second to the fifth energy levels of atomic oxygen. 

The electron-impact transition rates of the molecular species are compared to rates 
presented by Cherni[83], Telulet[84], and Gorelov[85]. Example transition rates for 
molecular nitrogen, nitric oxide, and the molecular nitrogen ion are shown in Figure 
4.5, Figure 4.6, and Figure 4.7, respectively. Again, good agreement is shown 
between the sampled rates and the reported rates. 
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Figure 4.5: Comparison of sampled and reported electronic energy level transition 
rates from the second to the third energy levels of molecular nitrogen. 
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Figure 4.6: Comparison of sampled and reported electronic energy level transition 
rates from the first to the fourth energy levels of nitric oxide. 
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Figure 4.7: Comparison of sampled and reported electronic energy level transition 
rates from the first to the third energy levels of the molecular nitrogen ion. 

4.3.2 Heavy Particle Impact Excitation 

Transitions due to heavy particle collisions are similar to electron-impact transitions, 
but they are significantly less important for most applications, as discussed by 
Johnston[86], where he shows that the process is negligible for atomic and molecular 
species in lunar-return shock layers. The comparisons shown in this section are, for 
the most part, as favorable as those in the previous section for electron-impact 
transition rates. Comparisons are made with values reported in References [85, 87- 
93]. For example, the transitions of the molecular nitrogen ion when colliding with 
molecular nitrogen, as shown in Figure 4.8, is generally within an order of magnitude 
of the reported values (note the 3-4 order of magnitude in the spread of data for the 
level 1-2 transition). There is also good agreement for collisions of molecular 
nitrogen and nitric oxide for most transitions, as exemplified in Figure 4.9 and Figure 
4.10, respectively. 
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Figure 4.8: Comparison of sampled and reported electronic energy level transition 
rates for the molecular nitrogen ion when colliding with molecular nitrogen. 
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Figure 4.9: Comparison of sampled and reported electronic energy level transition 
rates for molecular nitrogen when colliding with atomic nitrogen. 
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Figure 4.10: Comparison of sampled and reported electronic energy level transition 
rates for nitric oxide when colliding with molecular nitrogen. 

However, there seem to be certain combinations of heavy-impact partners where the 
current phenomenological model can’t accurately reproduce the measured/calculated 
transition rates. These combinations generally include transitions where the collision 
partner is molecular nitrogen, as shown in Figure 4.1 1, and the transition is from the 
lower levels (i.e., 1-2 or even some times 2-3). 
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Figure 4.11: Comparison of sampled and reported electronic energy level transition 
rates for the nitric oxide when colliding with molecular nitrogen. 

The current model also appears to have some difficulty in reproducing electronic 
energy level transitions in collisions of molecular oxygen, as shown in Figure 4. 12. 
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However, even though these comparisons have been less than ideal, they are not 
completely unreliable considering the four to five order of magnitude spread in the 
rates previously reported from the literature. 
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Figure 4.12: Comparison of sampled and reported electronic energy level transition 
rates for molecular oxygen when colliding with atomic oxygen. 

4.3.3 Equilibrium conditions 

In order for the electronic energy level procedures to be acceptable, they need to 
reproduce equilibrium conditions after the simulation has been allowed to run for a 
sufficient amount of time to ensure that detailed balance has been enforced. The first 
requirement is for the sampled translational and electronic temperatures to remain 
very near the input equilibrium temperature. The sampled translational, electron, and 
electronic temperatures for all of the species are presented in Table B.l at 10,000, 
20,000, and 30,000 K. All of the temperatures depart from the set equilibrium value 
by less than 0.5%. 

The resultant population distribution for the first seven energy levels of all atomic and 
molecular species (where applicable) for a range of temperatures is compared to the 
Boltzmann distribution in Table B.2 through Table B.3 1 . Again, the simulation 
represents the equilibrium state very well, with sampled departures from the 
Boltzmann distribution less than 1% for most levels where the number fraction in the 
state is above 1.0x10°. The larger departure from the equilibrium distribution is 
expected for levels with a very small population due to statistical noise. 


4.4 Associative Ionization Reactions 


Because of their low threshold energies, associative ionization reactions play an 
important role in determining the level of ionization and structure of the flow field at 
hypersonic reentry conditions. There exists very little experimental or computational 
data in the literature regarding the associative ionization reactions in air. The 
associative ionization reaction N + O NO + e~ has the lowest energy threshold, 
however the associative ionization reactions of N + N <-» iV, + e~ and 
O + O <-» 0\ + e~ are also very important mechanisms in the current application. The 
sampled DSMC reaction rates of endothermic and exothermic reactions of the N + O 
reaction are compared to measured rates from the literature in Figure 4.13 and Figure 
4.14, respectively. 
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Figure 4.13: Endothermic associative reaction rates for N + O — > NO 1 + e 
compared to rates from the literature. 
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Figure 4.14: Exothermic associative reaction rates for NO + + e —> N + 0 compared 
to rates from the literature. 

In these figures, the black lines are the analytical solutions from Equations (3.2.3) and 
(3.2.6), the black symbols are the sampled rates from the DSMC code, and the first of 
the values from literature is the curve used by LAURA. As discussed earlier, there is 
a wide spread in the reaction rates reported in the literature. Given the assumed 
measurement uncertainty, the Q-K reaction rates are in good agreement with the rates 
from the literature. 

4.5 Ionization Reactions 

Another important reaction mechanism for ionized, hypersonic reentry problems is 
the ionization of an atom or molecule by direct impact. There is very limited data in 
the literature for heavy particle impact ionization and these reactions are not generally 
included in reentry computations, so they will not be discussed in this work. There is, 
however, an appreciable amount of data related to electron impact ionization 
reactions, especially for atomic species. The electron impact ionization reactions 
generally considered in reentry computations are those for atomic nitrogen and 
oxygen. Reaction rate curves for the electron impact ionization reaction for atomic 
oxygen are presented in Figure 4.15. 
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Figure 4. 15: Endothermic ionization reaction rates for O + e —> O' + e + e 
compared to rates from the literature. 

4.6 Charge Exchange Reactions 

In charge exchange reactions, a singly charged ion (such as or O ) collides with a 
neutral atom or molecule and captures one of its electrons, thereby becoming a 
neutral atom. Only a small amount of energy is transferred to the electron donor, so 
the newly created neutral particle retains most of the original energy of the ion. 

These can therefore be an important mechanism for the creation of high-energy 
neutral particles. The measured reaction rates from the endothennic and exothermic 
charge exchange reactions corresponding to N 2 + N 1 —> N 2 + N are shown in Figure 
4.16 and Figure 4.17, respectively, compared to the values from the literature as 
shown. 
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Figure 4.16: Endothermic charge exchange reaction rates for N + + N 2 — » N 2 + N 
compared to rates from the literature. 
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Figure 4. 17: Exothermic charge exchange reaction rates for N 2 + N — > N + + N 2 
compared to rates from the literature. 

4.7 Exchange Reactions Involving Charged Species 

The remainder of relevant chemical reactions involves various exchange reactions 
involving charged species. As an example, the endothermic and exothermic reaction 
rates for the reaction NO + + N <-> N 2 + O are shown in Figure 4. 18 and Figure 4. 19, 
respectively, compared to values from the literature as shown. 
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Figure 4. 18: Endothermic reaction rate for the exchange reaction involving charged 
species NO + + N <-» + O with comparisons to rates from the literature. 
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Figure 4.19: Exothermic reaction rate for the exchange reaction involving charged 
species AC + O NO 1 + N with comparisons to rates from the literature. 
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Chapter 5: Modeling of the FIRE II Reentry Flow with Ionization and 

Electronic Energy Levels with DSMC 

In this chapter, the newly developed DSMC models presented in previous chapters 
are exercised on a demanding near-continuum reenty flow field with ionization. 

These results are compared with results from traditional DSMC models, a solution of 
the FIRE II flow field obtained using current best practices in computational fluid 
dynamics (CFD), and values inferred from flight. Both flow field and surface 
parameters are examined. To this end, three DSMC cases are summarized in Table 
5.1. 

Table 5.1: List of DSMC FIRE II cases. 


Case 

Electronic Energy 
Model 

Chemistry Model 

Case 1 

None 

TCE 

Case 2 

Present 

TCE 

Case 3 

Present 

Q-K 


5.1 The Project FIRE Flight Experiment 

The Project FIRE (Flight Investigation Reentry Environment) was an Apollo-era 
flight experiment to measure the radiative and convective heating on spacecraft 
material during atmospheric reentry at lunar return speeds (1 1.35 km/s)[21]. Project 
FIRE consisted of two sub-orbital flight experiments launched April 14, 1964 and 
May 22, 1965 from Cape Kennedy, Florida. For the current study, computations are 
performed for the reentry conditions of the second flight at 1634 seconds, which 
corresponds to an altitude of 76.42 km, and conditions for the simulation are provided 
in Table 5.2. The geometry and dimensions of the second vehicle in the experimental 
series, FIRE II, are presented in Figure 5.1. The forebody of the FIRE II reentry 
vehicle consisted of three phenolic-asbestos heat shields sandwiched between 
beryllium calorimeters. The first two heat shield and calorimeter packages were 
designed to be ejected after the onset of melting, yielding heating data free of the 
effects of ablation and three separate data gathering periods. Calorimeter plugs and 
radiometers were located at various positions on the heat shields. The flight condition 
examined in this study occurred during the first data collection period. 


Table 5.2: Free stream and surface parameters for FIRE II entry condition. 


Time (s) 

Altitude (kg/m 3 ) n x (m 3 ) 

(km) 

Foe (m/s) 

Te(K) 

Kn D Twaii (K) X N2 X 0 2 

1634 

76.42 3.72E-05 7.73e20 

11360 

195 

2.59E-03 615 0.76 0.24 
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Figure 5.1: Schematic of FIRE II test vehicle. 

5.2 Details of the Continuum Flow Model and Non-Equilibrium Radiation 
Model 

Two different numerical approaches are applied: a CFD solution of the Navier-Stokes 
equations using the Langley Aerothermodynamic Upwind Relaxation Algorithm code 
and particle -based DSMC computations using the DAC code, which is the focus of 
this study. All computations presented herein were steady state and axisymmetric 
with the x-direction along the axis. In each case, an 1 1-species, 21-reaction, 
thermochemical non-equilibrium approach is adopted. The surface catalytic 
boundary condition used in the CFD solution was “fully-catalytic,” which means that 
all atomic species recombine to their molecular counterpart (i.e., N + N — > N 2 ) and 
also charged heavy species recombine with electrons to fonn their charge neutral 
counterpart (i.e., N + + e~ — > N). 

The Langley Aerothennodynamic Upwind Relaxation Algorithm (LAURA) is a high 
fidelity, structured grid computer code, specialized for hypersonic re-entry physics, 
utilizing state-of-the-art algorithms for computational fluid dynamic simulations[94, 
95]. Key elements of LAURA include Roe’s averaging[96] and Yee’s Symmetric 
Total Variation Diminishing (STVD)[97] formulation of second-order, inviscid flux. 
Yee’s STVD fonnulation has been found to be exceptionally robust and Courant 
number independent using first point-implicit and then line-implicit relaxation for 
hypersonic flow simulations. 

A two-temperature thermochemical non-equilibrium model[95] is applied. The 
chemical reaction rates are compiled from previous studies (listed as the references 
for the first reaction rate parameters for each reaction in Table A. 14). The 
thermophysical properties are taken from the work of Mcbride et a/.[98] and Gupta et 





al. [76]. Multicomponent diffusion is approximated using Sutton and Gnoffo’s[99] 
approximate-corrected approach. 

The shock-layer radiation is modeled with the HARA (High-temperature 
Aerothermodynamic RAdiation) code. The details of this code for treating air species 
are presented by Johnson et al. [19, 20], Briefly, it is based on a set of atomic levels 
and lines obtained from the National Institute of Standards and Technology (NIST) 
online database[100] and the Opacity Project[101], as well as atomic bound-free cross 
sections from the TOPbase[102], The negative nitrogen and oxygen ions are treated 
using cross sections suggested by Soon and Kunc[103] and Chauveau et a/.[104], 
respectively. The molecular band systems are treated using a smeared-rotational band 
(SRB) model[105], which was shown by Johnston et al.[ 19] to be sufficient for 
treating VUV absorbing and optically-thin emitting band systems in air. The 
molecular data for modeling these band systems are obtained from Laux[106], except 
for the VUV N 2 systems, which are obtained from various other sources [107- 109]. 
The non-Boltzmann modeling of the atomic and molecular electronic states is based 
on a set of electron-impact excitation rates compiled from the literature and presented 
in detail by Johnston et al.[ 20], Following the work of Park[l 10], the quasi-steady 
state assumption is made when solving the Master Equation. The tangent-slab 
approximation is applied to calculate radiative flux and the divergence of the radiative 
flux. 


5.3 DSMC Simulation Parameters 

In this section, the parameters used in the simulation of the FIRE II flight experiment 
are discussed. 

5.3.1 Collision Modeling 

The baseline collision model implemented in DAC, for both neutral and charged 
species interactions, is that originally published by Bird[41], in which all particle 
interactions are treated using the VHS collision model. The reference diameters of 
the heavy ions O2 f , N} + , O , N + , and NO are set equal to the diameters of their 
associated neutral particles, and the reference diameter of the electrons is set to le-10 
m. This value of the electron reference diameter is much larger than the classical 
value of 6e-15 m, and is used to produce electron collision cross sections of the 
required magnitude within the mathematical framework of the VHS collision model. 
A complete list of the VHS reference diameters for all species and the parameters T re f 
and mis given in Table A.l. As shown previously, the present choice of electron 
reference diameter does not introduce significant error in the DSMC results. 

5.3.2 Internal Energy Modeling 

For the current thesis, the rotational collision number is set to unity. The reason for 
this is because the solutions are being compared to a computational fluid dynamics 
(CFD) solution where the translational and rotational temperatures are assumed equal. 
Therefore, for DSMC to best match CFD, the rotational energy is sampled from the 



equilibrium distribution at each collision (Z rot = 1.0). The vibrational energy 
relaxation model used is that based on the work of Millikan and White and the 
parameters used for each of the molecular species are presented in Table A. 2. 

One goal of this work is to examine the effects of adding an electronic energy 
transition algorithm to DSMC solutions of hypersonic (re)entry flows. To that end, 
the electronic energy levels of all relevant simulated species (O 2 , N 2 , O, N, NO, 0 2 + , 
N 2 + , 0 + , N + , and NO + ) are included in Cases 2 and 3. The energy levels included can 
be viewed in Table A. 4 through Table A. 13. The energy levels for O and N were 
taken from Johnston[86], while the remaining species were taken from Hash[l 1 1]. A 
more comprehensive set of energy levels could be compiled for 0 + and N + , but these 
are the energy levels used by the CFD code to be compared with, as well as the 
nonequilibrium radiation code. 

5.3.3 Chemistry Modeling 

The baseline chemistry model is the TCE model using the reaction rates used in 
LAURA. These rates are listed in Table A. 14 as the first rate listed for each reaction. 
As discussed in Section 2.2.3, the exothermic (reverse) reaction rates are then 
determined by considering the statistical mechanics equilibrium constant, and then 
the coefficients are chosen from those in the literature to most closely match the rate 
as computed by statistical mechanics. 

The other chemistry model, which is one of the main contributions, is the Q-K 
chemistry model and the extensions to charged species. The characteristic 
temperatures and reference parameters for vibrational relaxation can be found in 
Table A.l and Table A. 3, respectively. Additionally, the parameters required to 
account for the statistical mechanics equilibrium constant and adjust the activation 
energies for recombination and exchange/charge exchange can be found in Table 
A. 16 and Table A. 17, respectively. Comparisons of the Q-K reaction set to various 
rates found in the literature can be viewed in Appendix C. 

5.3.4 Electric Field Modeling 

The electric field model implemented in the current study, which was also added to 
DAC, is a hybrid of Bird’s model[26] and that of Boyd[30], In this model, the 
electron particles are moved along with their associated ion during the movement 
portion of the DSMC algorithm. However, the electrons retain their individual 
velocity components, which are used in the collision portion of the DSMC algorithm. 
Neither the electrons nor ions actually experience an electric field, in that their 
velocity components are not adjusted to account for an electric field as one is not 
computed. 

5.3.5 Gas-Surface Interactions 

The surface catalycity must be treated differently with the present DSMC simulations 
because the CFD solutions use the fully-catalytic surface boundary condition to 
account for the higher surface temperatures in the continuum regime. However, there 



is evidence[25] that, for surface temperatures such as those at the higher altitudes for 
the FIRE II experiment, the catalycity of neutral atomic species recombining is much 
closer to zero. Therefore, the interaction of particles with the vehicle surface is 
managed by assuming both a fully catalytic surface as was used with the CFD 
solution for comparison and also a partially catalytic surface that only allows charged 
species to become neutralized. The fully catalytic simulations will be referred to 
simply by their case number (i.e., Case 1) and the partially catalytic simulations will 
have an extra notation indicating that it is partially catalytic (i.e., Case l pc ). The 
remaining collisions (neutral molecules), along with the recombined particles, are 
then allowed to reflect diffusely from the surface, after thermally accommodating to 
the specified wall temperature. 

5.4 FIRE II Simulations 

In this section, a progression of DSMC solutions is discussed and compared with 
CFD beginning with a traditional DSMC solution without electronic energy levels 
(Case 1), followed by the addition of electronic energy levels (Case 2), concluding 
with the Q-K chemistry model (Case 3). At the end of the section, a discussion of the 
surface heating, both convective and radiative, is included where the DSMC rates are 
compared to both the CFD solution and values inferred from flight. 

5.4.1 Baseline Solution 

The baseline solution, referred to as Case 1 in Table 5.1, is the FIRE II simulation at 
1634 seconds using the TCE chemistry model without the addition of electronic 
energy levels. Various temperatures for this case are compared to that from the CFD 
solution in Figure 5.2. Clearly, the exclusion of electronic energy levels from the 
DSMC solution changes the shock-layer profile. The shock standoff distance of the 
DSMC solution is nearly 40% greater than that of the CFD solution, while the shock- 
layer temperature is also greater than that predicted by CFD. It can be seen that, 
inside the shock layer of the DSMC solution, the temperatures are not in equilibrium. 
This has been reported in Refs. [25] and [4] as an effect of the addition of ionization 
physics at this altitude. There is some scatter present in the data for both the 
rotational and vibrational temperatures, which may be caused by the low number 
density of molecular species in the flow field. The spike in vibrational temperature 
outside of the shock layer is the result of trace species of vibrationally excited 
molecules that have escaped from the shock layer. This is a phenomenon typical of 
flows such as this in DSMC simulations. Since the physics are clearly not the same, 
comparisons of species concentrations are not provided for this case. 
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Figure 5.2: Translational temperature comparison between CFD and Case 1 (no 
electronic energy levels). 

5.4.2 Addition of Electronic Energy Levels 


Electronic energy levels were included in DSMC Case 2 and are more appropriate for 
comparison to the CFD solution. The various temperature profiles from CFD and 
Case 2 are compared in Figure 5.3. The DSMC shock standoff distance has been 
reduced as compared to Case 1 and now the location and magnitude of the peak 
translational temperatures of the CFD and DSMC Case 2 solutions compare quite 
well. The differences between the solutions are the larger equilibrium temperature 
within the shock layer and the increased electronic temperature that persists upstream 
of the location of maximum translational temperature, which is discussed below. 
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Figure 5.3: Comparison of temperature profiles between the CFD solution and Case 
2 (electronic energy levels added). 

Comparisons of species mole fractions between the CFD solution and Case 2 are 
presented in Figure 5.4 and Figure 5.5 for neutral and charged species, respectively. 
In general, the species fractions compare well. As with the vibrational and rotational 
temperatures, there is noise in the mole fractions extracted from the DSMC 
simulation below around 0.05% because of the low sample size of these species. The 
DSMC solution contains trace species further upstream of the shock front as 
compared to the CFD solution. This is most likely caused by rarefaction effects 
where there are insufficient collisions to keep all particles within the shock layer. 

This also accounts for the increased electronic temperatures upstream of the 
maximum translational temperature as seen in Figure 5.3, and is analogous to what 
occurs with the vibrational temperature upstream of the shock, except with the 
electronic energy levels, not only do the molecules carry energy upstream of the 
shock, but also the atomic species. Even for a small population of excited particles, 
the electronic temperature can be large. 
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Figure 5.4: Comparison of neutral species mole fractions between the CFD solution 
and Case 2. 



Figure 5.5: Comparison of charged species mole fractions between the CFD solution 
and Case 2. 

5.4.3 Quantum-Kinetic Solution 

As was discussed earlier, DSMC Case 3 switches chemistry modules from the 
traditional TCE model to the new Q-K model. The various temperature profiles from 
CFD and Case 3 are compared in Figure 5.6. For this case, the peak translational 
temperature is slightly lower and moved upstream as compared to the CFD solution. 
This discrepancy may attributable both to the difference in the resultant equilibrium 
reaction rates discussed in Chapter 4 and also because the Q-K model does not rely 


1 1 


on equilibrium-based information. If the Q-K reaction rates are compared to the 
reaction rates used in the CFD solution (the first referenced rate in each figure) in 
Appendix C, there can be a large discrepancy, even though the Q-K rate is generally 
still well within the spread of reported reaction rates. 
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Figure 5.6: Comparison of temperature profiles between the CFD solution and Case 
3 (with electronic energy levels and Q-K chemistry). 


Comparisons of species mole fractions are presented between the CFD solution and 
DSMC Case 3 in Figure 5.7 and Figure 5.8 for neutral and charged species, 
respectively. Again, the mole fractions generally compare well between the CFD and 
Case 2 solutions with some notable differences. One difference is again the presence 
of trace species well upstream of the shock layer. 
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Figure 5.7: Comparison of neutral species mole fractions between the CFD solution 
and Case 3. 



Figure 5.8: Comparison of charged species mole fractions between the CFD solution 
and Case 3. 


5.4.4 Comparison of Surface Heating 

The FIRE II calorimeter data presented by Comette[21] provide heating values that 
contain the convective heating plus the contribution of the radiative flux absorbed by 
the beryllium calorimeter. In order to compare to these data, the code HARA was 
used to post-process the flow field solutions to obtain the non-equilibrium radiation 
heating to the surface. It was determined[21] that the fraction of the radiative flux 
absorbed by the calorimeter, a, was 0.72. 
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The computed heat flux to the surface of the vehicle is compared to the values 
inferred from the flight data at a radial location of approximately 0.02 m in Table 5.3. 
When comparing the convective heating component between the CFD and fully 
catalytic DSMC solutions including electronic energy levels, very good agreement is 
found. The solutions agree to within 1% of each other. However, DSMC Case 1, 
without electronic energy levels, has a convective heating rate approximately 1 5 
W/cm 2 higher, again highlighting the fact that electronic energy levels are required 
for the simulation. When the cases with electronic energy levels are compared to the 
inferred flight value, the convective heating is over-predicted by 34%. If the partially 
catalytic solutions are then considered, they under-predict the convective heating by 
around 15%. Therefore, the surface of the flight vehicle was likely at least partially 
catalytic to atomic recombination. 

The radiative heating component is also listed in Table 5.3 in the third column. Only 
the DSMC simulations with electronic energy levels modeled have been post- 
processed for radiative heating. Because of the differences in molar fractions 
between the DSMC simulations and the CFD solution, there is about a 10 W/cm 2 
increase in radiative heating for the DSMC simulations. There are generally more 
species present in the shock layer that radiate in the DSMC simulations, such as N, O, 
and N 2 + . The increase in radiative heating for the partially catalytic DSMC 
simulations is accounted for by the increased levels of N and O near the surface since 
these species do not recombine upon colliding with the wall. 

Table 5.3: FIRE II surface heating at r = 0.02 m. 
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Chapter 6: Conclusions 

This chapter contains a summary of the results of this thesis, the original 
contributions made to the field, and identifies future research directions. 

6,1 Summary of Results 

Extensions to the recently introduced kinetic-theory approach for computing chemical 
reaction rates (the Q-K method) have been developed and validated for DSMC for 
electronic energy level transitions and reactions involving charged species. These 
new models make use of the principles of microscopic reversibility and molecular- 
level energy exchange to predict the probability that an energy level transition or a 
chemical reaction occurs during a collision between two particles. Procedures have 
been defined that are required to implement these models in the DSMC technique. 

For each of these new Q-K reaction-rate models, comparisons have been made 
between the predicted analytic curves, the sampled transition or reaction rates from 
DSMC simulations, and rates obtained from the literature. The analytic and sampled 
values compare favorably to the rates from the literature; the historical rates often 
differ by several orders of magnitude for the same reactions. It should be noted 
again, however, that the Q-K model is phenomenological in nature and there may be 
cases for which it leads to rates that are outside of the range of uncertainty of the 
measured rates. Even so, the usefulness of these models has been shown as 
engineering models, especially where there are no data or the uncertainty in reported 
values is quite large. If there are reactions for which there is a high degree of 
certainty in the reported rate, they can be specified to use the TCE model. 

The proposed extensions to the Q-K model were implemented in the DSMC Analysis 
Code (DAC) of NASA and tested on the Project Fire II flight experiment. The 
addition of electronic energy levels to the standard chemistry model resulted in better 
agreement between DSMC and CFD. Both the location and magnitude of the 
maximum translational temperature agree. The DSMC solution without electronic 
energy levels had predicted a higher peak translational temperature and shock 
thickness approximately 40% greater than the CFD solution. However, because of 
rarefaction effects, trace species with elevated electronic energies were found 
upstream of the shock layer, resulting in an elevated electronic temperature. Similar 
results were found when the chemistry model was changed from the TCE model to 
the Q-K model. The only major difference was a slightly thicker shock layer and 
lower peak translational temperature because of the differences in resultant reaction 
rates. Most likely, the Q-K model better represents the reaction rates because it is a 
phenomenological approach, whereas for the other models, the reaction rates are 
extrapolated from low temperature experiments. 

Once surface convective heating had been obtained, the flow field was post-processed 
using the HARA non-equilibrium radiation solver to compare DSMC, CFD and 
values obtained from the flight experiment. The convective heating of the fully 
catalytic DSMC simulations with electronic energy levels included compared to 
within 1% of values obtained using CFD. These values are, however, around 34% 



higher than the value inferred from flight. When a partially catalytic surface 
boundary condition is employed, allowing only electron/ion recombinations, the 
convective heating dropped to levels on the order of 15% lower than the flight value. 
Therefore, the fully catalytic and partially catalytic simulations have been shown to 
bound the flight value of convective heating. The radiative heating computed from 
the DSMC simulations is around 10 W/cni 2 higher than that of the CFD simulation, 
which is caused by slightly larger amounts of N, O, and N 2 + . However, the 
nonequilibrium radiation solver is required to make certain assumptions in its 
calculations, so a direct input of the sampled non-Boltzmann distribution of electronic 
energy levels from the DSMC simulation would negate the necessity of the 
assumptions. 

6.2 Contributions 

The ability to simulate rarefied, ionized, hypersonic flows has been added to the suite 
of simulation capabilities for the National Aeronautics and Space Administration. 
Specifically, physical models for the treatment of charged particles using the DSMC 
method have been implemented in DAC. It has been shown that the inclusion of 
electronic energy levels, which few other DSMC solvers have the capability to do, is 
necessary for lunar-reentry velocity missions for the accurate prediction of shock 
shape and surface heating. The ability to predict electronic energy level transitions 
and resulting distributions using only energy level data has been developed as a result 
of this thesis. Also, the QK chemistry model was extended to include reactions 
involving charged species. This model allows the inclusion of chemical reactions for 
which specific reaction rate data is not known and may be more appropriate for non- 
equilibrium conditions. Both electronic energy level transition rates and chemical 
reaction rates have been shown to agree well with rates from the literature. 
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Appendix A: Species and Chemistry Data 


Table A. 1 : Parameters used in the VHS molecular model and characteristic 
temperatures. 



1 3 0.0000E+01 

2 2 1.5728E-19 a \ g 

3 1 2.621 IE-19 b l X + g 

4 1 6.5663E-19 C l X~ 

5 6 6.8912E-19 A' 3 A u 

7.0306E-19 A 3 Xl 

q ^ 


6 


3 












Level g E e i (J) State 

7 3 9.9267E-19 B 3 ^ 


Table A. 5: Electronic energy levels of N 2 . 


Level g E ei (J) State 


1 

1 

0.0000E+01 

x's; 

2 

3 

9.9727E-19 

a 3 e: 

3 

6 

1.1843E-18 

f? 3 n g 

4 

6 

1 . 1 88 IE- 1 8 

W 3 A„ 

5 

3 

1.3 165E-18 

b ,3 k 

6 

1 

1.3538E-18 

n'X 

7 

2 

1.3763E-18 


8 

2 

1.4483E-18 

w'A,, 

9 

5 

1.5415E-18 

A ' 5 K 

10 

1 

1.6925E-18 


11 

6 

1.7242E-18 


12 

6 

1.7707E-18 

C 3 n. 

13 

10 

1.8474E-18 

c" 5 n„ 

14 

6 

1.9388E-18 

c' 3 n„ 

15 

6 

2.0778E-18 

h 3 o„ 


Table A. 6: Electronic energy levels of O. 


Level 

g 

E(J) 

Core 

Config. 

Term 

1 

9 

1.55E-21 

2s 2 

2p 4 

3 P 

2 

5 

3.15E-19 

2s 2 

2p 4 

‘d 

3 

1 

6.71E-19 

2s 2 

2p 4 

'S 

4 

5 

1.47E-18 

2s 3 2p 2 

3s 

5 S 0 

5 

3 

1.53E-18 

2s 3 2p 2 

3s 

3 S° 

6 

15 

1.72E-18 

2s 

3p 

5p 

7 

9 

1.76E-18 

2s 3 2p 2 

3p 

3 P 

8 

5 

1.90E-18 

2s 3 2p 2 

4s 

5 S ° 

9 

3 

1.91E-18 

2s 3 2p 2 

4s 

3 S° 

10 

25 

1.94E-18 

2s 3 2p 2 

3d 

S D° 

11 

15 

1.94E-18 

2s 3 2p 2 

3d 

3 d° 

12 

15 

1.97E-18 

2s 3 2p 2 

4p 

5p 

13 

9 

1.98E-18 

2s 3 2p 2 

4p 

3 P 

14 

15 

2.01E-18 

2s 3 2p 2 

3s 

3 d° 

15 

5 

2.03E-18 

2s 3 2p 2 

5s 

5 S ° 

16 

3 

2.03E-18 

2s 3 2p 2 

5s 

3 S° 

17 

5 

2.04E-18 

2s 3 2p 2 

3s 

‘d 0 

18 

25 

2.04E-18 

2s 3 2p 2 

4d 

5 d° 









Level 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

g 

15 

35 

21 

15 

9 

25 

15 

35 

21 

288 

392 

512 

648 

800 

E(J) 

2.04E-18 

2.05E-18 

2.05E-18 

2.06E-18 

2.06E-18 

2.09E-18 

2.09E-18 

2.09E-18 

2.09E-18 

2.12E-18 

2.14E-18 

2.15E-18 

2.15E-18 

2.16E-18 

Core 

2s 3 2p 2 

2s 3 2p 2 

2s 3 2p 2 

2s 3 2p 2 

2s 3 2p 2 

2s 3 2p 2 

2s 3 2p 2 

2s 3 2p 2 

2s 3 2p 2 

2s 3 2p 2 

2s 3 2p 2 

2s 3 2p 2 

2s 3 2p 2 

2s 3 2p 2 

Config. 

4d 

4f 

4f 

5p 

5p 

5d 

5d 

5f 

5f 

n=6 

n=7 

n=8 

n=9 

n=10 

Term 

3 D° 

5 f 

3 f 

5p 

3 P 

5 d° 

3 d° 

5 f 

3 f 

Table A. 7: Electronic energy levels of N. 

Level 

g 

E d { J) 

Core 

Config. 

Term 

1 

4 

0.0000E+01 

2s 2 

2p 3 

4 S° 

2 

10 

3.8 1 95E-19 

2s 2 

2p 3 

2 d° 

3 

6 

5.7287E-19 

2s 2 

2p 3 

2pO 

4 

12 

1.6554E-18 

2s 2 2p 2 

3s 

4 p 

5 

6 

1.7122E-18 

2s 2 2p 2 

3s 

2 p 

6 

12 

1.7507E-18 

2s 

2p 4 

4 p 

7 

2 

1.8589E-18 

2s 2 2p 2 

3p 

2 S° 

8 

20 

1.8839E-18 

2s 2 2p 2 

3p 

4 d° 

9 

12 

1.8973E-18 

2s 2 2p 2 

3p 

4po 

10 

4 

1.9219E-18 

2s 2 2p 2 

3p 

4 S° 

11 

10 

1.9235E-18 

2s 2 2p 2 

3p 

2 d° 

12 

6 

1.9426E-18 

2s 2 2p 2 

3p 

2po 

13 

10 

1.9798E-18 

2s 2 2p 2 

3s 

2 d 

14 

12 

2.0598E-18 

2s 2 2p 2 

4s 

4 P 

15 

6 

2.0698E-18 

2s 2 2p 2 

4s 

2 P 

16 

6 

2.0784E-18 

2s 2 2p 2 

3d 

2 P 

17 

28 

2.0802E-18 

2s 2 2p 2 

3d 

4 f 

18 

14 

2.0828E-18 

2s 2 2p 2 

3d 

2 f 

19 

12 

2.0828E-18 

2s 2 2p 2 

3d 

4 P 

20 

20 

2.0859E-18 

2s 2 2p 2 

3d 

4 d 

21 

10 

2.0884E-18 

2s 2 2p 2 

3d 

2 d 

22 

2 

2.1 151E-18 

2s 2 2p 2 

4p 

2 S° 

23 

20 

2.1220E-18 

2s 2 2p 2 

4p 

4 d° 

24 

12 

2.1258E-18 

2s 2 2p 2 

4p 

4po 

25 

10 

2.1300E-18 

2s 2 2p 2 

4p 

2 d° 

26 

4 

2.1343E-18 

2s 2 2p 2 

4p 

4 S° 

27 

6 

2.1377E-18 

2s 2 2p 2 

4p 

2po 

28 

90 

2.191 2E- 1 8 

2s 2 2p 2 

4d i P,4d 4 F,4d 4 D,4d 4 F,4d 4 P,4d~D 

29 

126 

2.1946E-18 

2s 2 2p 2 

4f 4 D,4f 4 F,4f 4 G,4f 2 D,4f ! F,4f 2 G 

30 

450 

2.2368E-18 

2s 2 2p 2 

n= 

: 5 

31 

648 

2.2703E-18 

2s 2 2p 2 

n= 

6 

32 

822 

2.2864E-18 

2s 2 2p 2 

n= 

■1 

33 

1152 

2.2968E-18 

2s 2 2p 2 

n= 

: 8 

34 

1458 

2.3040E-18 

2s 2 2p 2 

n= 

9 






Level g E et (J) Core Config. Term 

35 1800 2.3091E-18 2s z 2p 2 n=10 


Table A. 8: Electronic energy levels of NO. 


Level 

g 

E e l( J) 

State 

1 

4 

0.0000E+01 

x 2 n r 

2 

8 

7.5485E-19 

a 4 n, 

3 

2 

8.7218E-19 

a 2 e + 

4 

4 

9.1 1 15E-19 

b 2 n r 

5 

4 

9.5349E-19 

b 4 YF 

6 

4 

9.7336E-19 

z/ 2 ® 

7 

4 

1.0343E-18 

c 2 n 

8 

2 

1.0533E-18 

d 2 f 

9 

4 

1.1979E-18 

B' 2 A 

10 

2 

1.2032E-18 

e 2 x + 

11 

4 

1.2269E-18 

F 2 A 

12 

4 

1.2401E-18 

// 2 n 

13 

2 

1.2410E-18 

H' 2 X + 

14 

2 

1.2485E-18 

g 2 z 

15 

2 

1 .25 15E-18 

/ 2 Z + 

16 

4 

1.2580E-18 

L 2 n 


Table A. 9: Electronic energy levels of 02 + . 


Level g E el (J) 

State 


1 

4 

0.0000E+01 

X 2 n g 

2 

8 

6.5380E-19 

« 4 n„ 

3 

4 

8.0594E-19 

A 2 n„ 

4 

6 

8.0650E-19 

6 ^ + 
u 

5 

4 

8.6013E-19 


6 

2 

9.2966E-19 

B 2 K 

7 

4 

9.8329E-19 


8 

4 

1.0568E-18 

c 2 o„ 

9 

4 

1.2177E-18 

d 2 n„ 

10 

4 

1.2276E-18 


11 

8 

1.301 IE-18 


12 

4 

1.311 IE- 1 8 

F 2 n g 

13 

2 

1.3243E-18 

G 2 *- g 


QQ 









Level 

8 

Eel (J) 

State 

14 

2 

1.3786E-18 

H 2 K 

15 

4 

1.4302E-18 

e A K 


Table A. 10: Electronic energy levels of N 2 + . 


Level g E e i (J) State 


1 

2 

0.0000E+01 

a 2 e; 

2 

4 

1.821 IE-19 

A 2 n„ 

3 

2 

5.0578E-19 

b 2 K 

4 

4 

5.0654E-19 

a A K 

5 

8 

8.2636E-19 

b4U s 

6 

8 

9.1377E-19 

c 4 n„ 

7 

4 

1.0492E-18 

c 2 n g 

8 

4 

1.0528E-18 

d 4 z;, 

9 

4 

1.1323E-18 

4 v + 

10 

4 

1 . 1 62 IE- 1 8 

D 2 O,, 

11 

8 

1.1919E-18 

f 4 n„ 

12 

8 

1.2316E-18 

S 4 A„ 

13 

4 

1.2713E-18 

E 2 \, 

14 

4 

1.2733E-18 

E 2 A„ 

15 

2 

1.2831E-18 

E 2 E U 

16 

2 

1.3309E-18 

E 2 A„ 

17 

4 

1.4302E-18 

E 2 A„ 


Table A. 1 1 : Electronic energy levels of 0 + . 


Level 

g 

Srf(J) 

Core 

Config. 

Term 

1 

4 

0.0000E+01 

2s 2 

2p 3 

4 S° 

2 

10 

5.3270E-19 

2s 2 

2p 3 

2 d° 

3 

6 

8.0386E-19 

2s 2 

2p 3 

2pO 


Table A. 12: Electronic energy levels ofN + . 


Level g E e i( J) Core Config. Term 


l 

9 

0.0000E+01 

2s 2 

2p 2 

3 P 

2 

5 

3.0424E-19 

2s 2 

2p 2 

'd 

3 

1 

6.4929E-19 

2s 2 

2p 2 

'S 

4 

15 

1.8324E-18 

2s 

2p 3 

3 d° 

5 

9 

2.1696E-18 

2s 

2p 3 

3po 

6 

5 

2.8642E-18 

2s 

2p 3 

‘d° 


QO 












Level g E el (J) Core Config. Term 

7 12 2.9610E-18 2s 2 2p 3s 3 P° 


Table A. 13: Electronic energy levels of NO + . 


Level 

g 

E el ( J) 

State 

1 

1 

0.0000E+01 

A‘E + 

2 

3 

1.0367E-18 

a 3 E + 

3 

6 

1 . 177 IE- 1 8 

b 3 n 

4 

6 

1.2293E-18 

c 3 A 

5 

3 

1.3457E-18 

d 3 YT 

6 

1 

1.381 4E- 1 8 

s'z 

7 

2 

1.4194E-18 

C 1 A 

8 

2 

1.4595E-18 

zyn 


Table A. 14: Endothermic chemical reaction rate coefficients (m 3 /molecule/s) used in 
the TCE model (the first rate listed for each reaction is the rate used in the 
simulations), kf = aT h cxp(-EJkT) 


Number 

Reaction 

a 

b 

C 

E a { J) 

E rx (J) 

Reference 

1 

O 2 + O 2 - 

■> 0 + 0 + O 2 

3.321 IE-09 

-1.50 

1.5 

8.1995E-19 

-8.1995E-19 

Park[63] 




4.5800E-1 1 

-1.00 


8.1970E-19 

-8.1970E-19 

Bird[ 15] 




5.9780E-12 

-1.00 


8.2150E-19 

-8.2150E-19 

Kang[64] 




5.3930E-1 1 

-1.00 


8.1970E-19 

-8.1970E-19 

Gimelshein[65] 

2 

O 2 + N 2 - 

-> 0 + 0 + N 2 

3.3211E-09 

-1.50 

1.5 

8.1995E-19 

-8.1995E-19 

Park[63] 




4.5800E-1 1 

-1.00 


8.1970E-19 

-8.1970E-19 

B ird [15] 




5.9780E-12 

-1.00 


8.2150E-19 

-8.2150E-19 

Kang[64] 




1.1980E-11 

-1.00 


8.1970E-19 

-8.1970E-19 

Gimelshein[65] 

3 

O 2 + 0 - 

->0 + 0 + 0 

1.6605E-08 

-1.50 

1.0 

8.1995E-19 

-8.1995E-19 

Park[63] 




1.3750E-10 

-1.00 


8.1970E-19 

-8.1970E-19 

B ird [15] 




5.9780E-12 

-1.00 


8.2150E-19 

-8.2150E-19 

Kang[64] 




3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa[9] 




1.4980E-10 

-1.00 


8.1970E-19 

-8.1970E-19 

Gimelshein[65] 




4.5700E-1 1 

-1.00 


8.1970E-19 

-8.1970E-19 

Johnston[66] 

4 

O 2 + N- 

-^0 + 0+N 

1.6605E-08 

-1.50 

1.0 

8.1995E-19 

-8.1995E-19 

Park[63] 




1.3750E-10 

-1.00 


8.1970E-19 

-8.1970E-19 

B ird [15] 




5.9780E-12 

-1.00 


8.2150E-19 

-8.2150E-19 

Kang[64] 




3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa[9] 




5.9930E-12 

-1.00 


8.1970E-19 

-8.1970E-19 

Gimelshein[65] 

5 

0 2 + NO - 

-> 0 + 0 + NO 

3.321 IE-09 

-1.50 

1.0 

8.1995E-19 

-8.1995E-19 

Park[63] 




4.5800E-1 1 

-1.00 


8.1970E-19 

-8.1970E-19 

B ird [15] 




5.9780E-12 

-1.00 


8.2150E-19 

-8.2150E-19 

Kang[64] 




5.9930E-12 

-1.00 


8.1970E-19 

-8.1970E-19 

Gimelshein[65] 

6 

O 2 + 02 + “ 

-> 0 + 0 + O 2 

3.321 IE-09 

-1.50 

1.5 

8.1995E-19 

-8.1995E-19 

Park[63] 




4.5800E-1 1 

-1.00 


8.1970E-19 

-8.1970E-19 

B ird [15] 




5.9780E-12 

-1.00 


8.2150E-19 

-8.2150E-19 

Kang[64] 




5.3930E-1 1 

-1.00 


8.1970E-19 

-8.1970E-19 

Gimelshein[65] 


on 








Number 

Reaction 

a 

b 

Q 

E a ( J) 

E„ (J) 

Reference 

7 

O2 + N2 + - 

■> 0 + 0 + N2 + 

3.321 IE-09 

-1.50 

1.5 

8.1995E-19 

-8.1995E-19 

Park[63] 




4.5800E-1 1 

-1.00 


8.1970E-19 

-8.1970E-19 

B ird [15] 




5.9780E-12 

-1.00 


8.2150E-19 

-8.2150E-19 

Kang[64] 




1.1980E-11 

-1.00 


8.1970E-19 

-8.1970E-19 

Gimelshein[65] 

8 

0 2 +0 + - 

■* 0 + 0 + 0 + 

1.6605E-08 

-1.50 

1.0 

8.1995E-19 

-8.1995E-19 

Park[63] 




1.3750E-10 

-1.00 


8.1970E-19 

-8.1970E-19 

B ird [15] 




5.9780E-12 

-1.00 


8.2150E-19 

-8.2150E-19 

Kang[64] 




3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa[9] 




1.4980E-10 

-1.00 


8.1970E-19 

-8.1970E-19 

Gimelshein[65] 




4.5700E-1 1 

-1.00 


8.1970E-19 

-8.1970E-19 

Johnston[66] 

9 

O2 + N + — 

-> 0 + 0 + N + 

1.6605E-08 

-1.50 

1.0 

8.1995E-19 

-8.1995E-19 

Park[63] 




1.3750E-10 

-1.00 


8.1970E-19 

-8.1970E-19 

Bird[ 15] 




5.9780E-12 

-1.00 


8.2150E-19 

-8.2150E-19 

Kang[64] 




3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa[9] 




5.9930E-12 

-1.00 


8.1970E-19 

-8.1970E-19 

Gimelshein[65] 

10 

O2 + NO - 

-+ 0 + 0 + NO + 

3.321 IE-09 

-1.50 

1.0 

8.1995E-19 

-8.1995E-19 

Park[63] 




4.5800E-1 1 

-1.00 


8.1970E-19 

-8.1970E-19 

Bird[ 15] 




5.9780E-12 

-1.00 


8.2150E-19 

-8.2150E-19 

Kang[64] 




5.9930E-12 

-1.00 


8.1970E-19 

-8.1970E-19 

Gimelshein[65] 

11 

N2 + O2 - 

-> N + N + O2 

1.1620E-08 

-1.60 

0.5 

1.5630E-18 

-1.5630E-18 

Park[63] 




6.1700E-09 

-1.60 


1.5610E-18 

-1.5610E-18 

B ird [15] 




3.1550E-13 

-0.50 


1.5600E-18 

-1.5600E-18 

Kang[64] 




4.9800E-08 

-1.60 


1.5630E-18 

-1.5630E-18 

Ozawa[9] 




7.9700E-13 

-0.50 


1 .561 0E- 1 8 

-1.561 0E- 1 8 

Boyd[67] 




3.1870E-13 

-0.50 


1 .561 0E- 1 8 

-1.5610E-18 

Gimelshein[65] 

12 

N2 -t - N2 — 

■> N + N + N 2 

1.1620E-08 

-1.60 

1.0 

1.5630E-18 

-1.5630E-18 

Park[63] 




6.1700E-09 

-1.60 


1 .561 0E- 1 8 

-1.5610E-18 

Bird[ 15] 




3.1550E-13 

-0.50 


1.5600E-18 

-1.5600E-18 

Kang[64] 




1.1600E-08 

-1.60 


1.5630E-18 

-1.5630E-18 

Ozawa[9] 




7.9700E-13 

-0.50 


1 .561 0E- 1 8 

-1.561 0E- 1 8 

Boyd[67] 

13 

N 2 + 0 - 

-^N + N + O 

4.9980E-08 

-1.60 

0.5 

1.5630E-18 

-1.5630E-18 

Park[63] 




1.8500E-08 

-1.60 


1 .561 0E- 1 8 

-1.5610E-18 

B ird [15] 




3.1550E-13 

-0.50 


1.5600E-18 

-1.5600E-18 

Kang[64] 




7.1400E-08 

-1.50 


1.5630E-18 

-1.5630E-18 

Boyd[67] 




3.1870E-13 

-0.50 


1 .561 0E- 1 8 

-1.5610E-18 

Gimelshein[65] 

14 

N 2 + N - 

-^N + N + N 

4.9980E-08 

-1.60 

1.0 

1.5630E-18 

-1.5630E-18 

Park[63] 




1.8500E-08 

-1.60 


1 .561 0E- 1 8 

-1.5610E-18 

B ird [15] 




6.7830E-08 

-1.50 


1.5600E-18 

-1.5600E-18 

Kang[64] 




7.1400E-08 

-1.50 


1.5630E-18 

-1.5630E-18 

Boyd[67] 




6.9000E-08 

-1.50 


1 .561 0E- 1 8 

-1.5610E-18 

Gimelshein[65] 

15 

N 2 + NO - 

-» N + N + NO 

1.1620E-08 

-1.60 

0.5 

1.5630E-18 

-1.5630E-18 

Park[63] 




6.1700E-09 

-1.60 


1.5610E-18 

-1.561 0E- 1 8 

Bird[ 15] 




3.1550E-13 

-0.50 


1.5600E-18 

-1.5600E-18 

Kang[64] 




4.9800E-08 

-1.60 


1.5630E-18 

-1.5630E-18 

Ozawa[9] 




7.9700E-13 

-0.50 


1 .561 0E- 1 8 

-1.561 0E- 1 8 

Boyd[67] 




3.1870E-13 

-0.50 


1.5610E-18 

-1.561 0E- 1 8 

Gimelshein[65] 

16 

N2 + 02 + - 

■> N + N + 02 + 

1.1620E-08 

-1.60 

0.5 

1.5630E-18 

-1.5630E-18 

Park[63] 




6.1700E-09 

-1.60 


1 .561 0E- 1 8 

-1.5610E-18 

B ird [15] 


Ol 



Number 

Reaction 

a 

b 


E a ( J) 

E rx (J) 

Reference 




3.1550E-13 

-0.50 


1.5600E-18 

-1.5600E-18 

Kang[64] 




4.9800E-08 

-1.60 


1.5630E-18 

-1.5630E-18 

Ozawa[9] 




7.9700E-13 

-0.50 


1 .561 0E- 1 8 

-1.5610E-18 

Boyd[67] 




3.1870E-13 

-0.50 


1 .561 0E- 1 8 

-1.5610E-18 

Gimelshein[65] 

17 

N 2 + N2 + — 

■> N + N + N 2 + 

1.1620E-08 

-1.60 

1.0 

1.5630E-18 

-1.5630E-18 

Park[63] 




6.1700E-09 

-1.60 


1 .561 0E- 1 8 

-1.5610E-18 

B ird [15] 




3.1550E-13 

-0.50 


1.5600E-18 

-1.5600E-18 

Kang[64] 




1.1600E-08 

-1.60 


1.5630E-18 

-1.5630E-18 

Ozawa[9] 




7.9700E-13 

-0.50 


1.5610E-18 

-1.561 0E- 1 8 

Boyd[67] 

18 

N 2 + o' — 

■> N + N + 0 + 

4.9980E-08 

-1.60 

0.5 

1.5630E-18 

-1.5630E-18 

Park[63] 




1.8500E-08 

-1.60 


1 .561 0E- 1 8 

-1.5610E-18 

B ird [15] 




3.1550E-13 

-0.50 


1.5600E-18 

-1.5600E-18 

Kang[64] 




7.1400E-08 

-1.50 


1.5630E-18 

-1.5630E-18 

Boyd[67] 




3.1870E-13 

-0.50 


1 .561 0E- 1 8 

-1.5610E-18 

Gimelshein[65] 

19 

N 2 + N + - 

■> N + N + N + 

4.9980E-08 

-1.60 

1.0 

1.5630E-18 

-1.5630E-18 

Park[63] 




1.8500E-08 

-1.60 


1 .561 0E- 1 8 

-1.5610E-18 

Bird[ 15] 




6.7830E-08 

-1.50 


1.5600E-18 

-1.5600E-18 

Kang[64] 




7.1400E-08 

-1.50 


1.5630E-18 

-1.5630E-18 

Boyd[67] 




6.9000E-08 

-1.50 


1 .561 0E- 1 8 

-1.5610E-18 

Gimelshein[65] 

20 

N 2 + NO + - 

■> N + N + NO + 

1.1620E-08 

-1.60 

0.5 

1.5630E-18 

-1.5630E-18 

Park[63] 




6.1700E-09 

-1.60 


1 .561 0E- 1 8 

-1.561 0E- 1 8 

B ird [15] 




3.1550E-13 

-0.50 


1.5600E-18 

-1.5600E-18 

Kang[64] 




4.9800E-08 

-1.60 


1.5630E-18 

-1.5630E-18 

Ozawa[9] 




7.9700E-13 

-0.50 


1 .561 0E- 1 8 

-1.561 0E- 1 8 

Boyd[67] 




3.1870E-13 

-0.50 


1.5610E-18 

-1.5610E-18 

Gimelshein[65] 

21 

N 2 + e — 

■> N + N + e' 

9.9670E-27 

2.60 

0.5 

1.5630E-18 

-1.5630E-18 

Bourdon[68] 




4.9800E-06 

-1.60 


1 .561 0E- 1 8 

-1.5630E-18 

Ozawa[9] 




1.4900E-05 

-1.60 


1 .561 0E- 1 8 

-1.561 0E- 1 8 

Park[69] 




1.3450E-10 

-1.28 


1.5630E-18 

-1.5630E-18 

Losev[70] 

22 

NO + 0 2 - 

-> N + 0 + O 2 

3.8300E-13 

-0.50 

1.0 

1.0420E-18 

-1.0420E-18 

Fujita[71] 




6.4760E-10 

-1.50 


1.0420E-18 

-1.0420E-18 

Kang[64] 




3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa[9] 




8.3000E-15 

0.00 


1.0420E-18 

-1.0420E-18 

Park[69] 




6.5900E-10 

-1.50 


1.0430E-18 

-1.0430E-18 

Gimelshein[65] 

23 

NO +N 2 - 

-> N + 0 + N 2 

3.8300E-13 

-0.50 

1.0 

1.0420E-18 

-1.0420E-18 

Fujita[71] 




6.4760E-10 

-1.50 


1.0420E-18 

-1.0420E-18 

Kang[64] 




3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa[9] 




8.3000E-15 

0.00 


1.0420E-18 

-1.0420E-18 

Park[69] 




6.5900E-10 

-1.50 


1.0430E-18 

-1.0430E-18 

Gimelshein[65] 

24 

NO + 0 - 

N + 0 + 0 

7.6600E-13 

-0.50 

1.0 

1.0420E-18 

-1.0420E-18 

Fujita[71] 




6.4760E-10 

-1.50 


1.0420E-18 

-1.0420E-18 

Kang[64] 




8.3000E-15 

0.00 


1.0420E-18 

-1.0420E-18 

Ozawa[9] 




1.8300E-13 

0.00 


1.0420E-18 

-1.0420E-18 

Park[69] 




1.3180E-08 

-1.50 


1.0430E-18 

-1.0430E-18 

Gimelshein[65] 

25 

NO + N- 

-> N + 0 + N 

7.6600E-13 

-0.50 

1.0 

1.0420E-18 

-1.0420E-18 

Fujita[71] 




6.4760E-10 

-1.50 


1.0420E-18 

-1.0420E-18 

Kang[64] 




3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa[9] 




1.8300E-13 

0.00 


1.0420E-18 

-1.0420E-18 

Park[69] 


cn 



Number 

Reaction 

a 

b 

Q 

E a ( J) 

E rx (J) 

Reference 




1.3180E-08 

-1.50 


1.0430E-18 

-1.0430E-18 

Gimelshein[65] 

26 

NO + NO - 

->■ N + 0 + NO 

3.8300E-13 

-0.50 

1.0 

1.0420E-18 

-1.0420E-18 

Fujita[71] 




6.4760E-10 

-1.50 


1.0420E-18 

-1.0420E-18 

Kang[64] 




3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa[9] 




8.3000E-15 

0.00 


1.0420E-18 

-1.0420E-18 

Park[69] 




1.3180E-08 

-1.50 


1.0430E-18 

-1.0430E-18 

Gimelshein[65] 

27 

NO + 02 + - 

— > N + 0 + 0 2 + 

3.8300E-13 

-0.50 

1.0 

1.0420E-18 

-1.0420E-18 

Fujita[71] 




3.8300E-13 

-0.50 


1.0430E-18 

-1.0430E-18 

B ird [15] 




6.4760E-10 

-1.50 


1.0420E-18 

-1.0420E-18 

Kang[64] 




3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa[9] 




8.3000E-15 

0.00 


1.0420E-18 

-1.0420E-18 

Park[69] 




6.5900E-10 

-1.50 


1.0430E-18 

-1.0430E-18 

Gimelshein[65] 

28 

NO + N 2 + - 

— > N + 0 + N 2 + 

3.8300E-13 

-0.50 

1.0 

1.0420E-18 

-1.0420E-18 

Fujita[71] 




3.8300E-13 

-0.50 


1.0430E-18 

-1.0430E-18 

Bird[ 15] 




6.4760E-10 

-1.50 


1.0420E-18 

-1.0420E-18 

Kang[64] 




3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa[9] 




8.3000E-15 

0.00 


1.0420E-18 

-1.0420E-18 

Park[69] 




6.5900E-10 

-1.50 


1.0430E-18 

-1.0430E-18 

Gimelshein[65] 

29 

NO + 0 + - 

-> N + 0 + 0 + 

7.6600E-13 

-0.50 

1.0 

1.0420E-18 

-1.0420E-18 

Fujita[71] 




7.6600E-13 

-0.50 


1.0430E-18 

-1.0430E-18 

Bird[ 15] 




6.4760E-10 

-1.50 


1.0420E-18 

-1.0420E-18 

Kang[64] 




8.3000E-15 

0.00 


1.0420E-18 

-1.0420E-18 

Ozawa[9] 




1.8300E-13 

0.00 


1.0420E-18 

-1.0420E-18 

Park[69] 




1.3180E-08 

-1.50 


1.0430E-18 

-1.0430E-18 

Gimelshein[65] 

30 

NO + N + - 

-> N + 0 + N + 

7.6600E-13 

-0.50 

1.0 

1.0420E-18 

-1.0420E-18 

Fujita[71] 




7.6600E-13 

-0.50 


1.0430E-18 

-1.0430E-18 

B ird [15] 




6.4760E-10 

-1.50 


1.0420E-18 

-1.0420E-18 

Kang[64] 




3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa[9] 




1.8300E-13 

0.00 


1.0420E-18 

-1.0420E-18 

Park[69] 




1.3180E-08 

-1.50 


1.0430E-18 

-1.0430E-18 

Gimelshein[65] 

31 

NO + NO + - 

-> N + 0 + NO + 

3.8300E-13 

-0.50 

1.0 

1.0420E-18 

-1.0420E-18 

Fujita[71] 




3.8300E-13 

-0.50 


1.0430E-18 

-1.0430E-18 

Bird[ 15] 




6.4760E-10 

-1.50 


1.0420E-18 

-1.0420E-18 

Kang[64] 




3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa[9] 




8.3000E-15 

0.00 


1.0420E-18 

-1.0420E-18 

Park[69] 




1.3180E-08 

-1.50 


1.0430E-18 

-1.0430E-18 

Gimelshein[65] 

32 

N + e'^ 

► N + + e" + e" 

5.8100E-08 

-1.00 

0.0 

2.3220E-18 

-2.3220E-18 

Park[63] 




4.1500E+03 

-3.82 


2.3300E-18 

-2.3300E-18 

Park[60] 




9.0000E-1 1 

-1.00 


1.5000E-18 

-1.5000E-18 

Ozawa[7] 




2.9890E-17 

0.60 


2.3300E-18 

-2.3300E-18 

Losev[70] 




5.8100E-08 

-1.00 


2.3300E-18 

-2.3300E-18 

Boyd[67] 




1.0000E-14 

0.00 


2.3300E-18 

-2.3300E-18 

Bird[ 15] 




1.8270E+02 

-3.14 


2.3330E-18 

-2.3300E-18 

Kang[64] 




8.4340E-14 

0.00 


2.3300E-18 

-2.3300E-18 

Farbar[5] 

33 

0 + e' -h 

► O 1 + e" + e" 

1.5900E-08 

-1.00 

0.0 

2.1880E-18 

-2.1880E-18 

Park[63] 




6.4800E+03 

-3.78 


2.1900E-18 

-2.1900E-18 

Park[60] 




3.0000E-10 

-1.00 


1.4000E-18 

-1.4000E-18 

Ozawa[7] 




8.6350E-18 

0.68 


2.181 0E- 1 8 

-2.1810E-18 

Losev[70] 


cn 



Number 

Reaction 

a 

b 

c 

E a ( J) 

E„ (J) 

Reference 




1.5900E-08 

-1.00 


2.181 0E- 1 8 

-2. 181 0E- 1 8 

Boyd[67] 




3.0000E-12 

0.00 


2.181 0E- 1 8 

-2.1810E-18 

Bird[ 15] 




5.9780E+01 

-2.91 


2.181 0E- 1 8 

-2.1810E-18 

Kang[64] 




1.0540E-14 

0.00 


2.181 0E- 1 8 

-2.1810E-18 

Farbar[5] 

34 

N 2 + 0 - 

-> NO + N 

9.9630E-17 

0.10 

0.0 

5.2470E-19 

-5.2470E-19 

Fujita[71] 




5.3000E-17 

0.10 


5.1750E-19 

-5.1750E-19 

Bird[ 1 5] 




1.1620E-16 

0.00 


5.2460E-19 

-5.2460E-19 

Kang[64] 




9.4500E-18 

0.42 


5.9280E-19 

-5.9280E-19 

Ozawa[9] 




1.0600E-12 

-1.00 


5.1750E-19 

-5.1750E-19 

Park[69] 




1.1200E-16 

0.00 


5.1750E-19 

-5.1750E-19 

Gimelshein[65] 

35 

NO + 0 - 

— » 0 2 + N 

1.3950E-17 

0.00 

0.0 

2.6790E-19 

-2.6790E-19 

Park[63] 




3.6000E-22 

1.29 


2.2380E-19 

-2.2380E-19 

Bird[ 1 5] 




5.3140E-21 

1.00 


2.7200E-19 

-2.7200E-19 

Kang[64] 




5.2790E-21 

1.00 


2.7190E-19 

-2.7190E-19 

Gimelshein[65] 

36 

N + N - 

* N 2 + + e 

7.3060E-23 

1.50 

0.0 

9.3190E-19 

-9.3190E-19 

Park[63] 




3.3870E-17 

0.00 


9.3200E-19 

-9.3200E-19 

Park[60] 




2.9800E-20 

0.77 


9.3400E-19 

-9.3400E-19 

Bird[ 1 5] 




2.3250E-17 

0.00 


9.3610E-19 

-9.3610E-19 

Kang[64] 

37 

N + O -> 

■ N0 + + e" 

8.8010E-18 

0.00 

0.0 

4.4040E-19 

-4.4040E-19 

Park[63] 




8.7660E-18 

0.00 


4.4000E-19 

-4.4000E-19 

Park[60] 




1.4600E-21 

1.00 


4.4000E-19 

-4.4000E-19 

Park[69] 




2.5500E-20 

0.37 


4.4220E-19 

-4.4220E-19 

Bird[ 1 5] 




2.3250E-24 

1.50 


4.4040E-19 

-4.4040E-19 

Kang[64] 




1.4990E-20 

0.50 


4.4730E-19 

-4.4730E-19 

Gupta[72] 

38 

0 + 0 - 

> 0 2 + + e 

1.1790E-27 

2.70 

0.0 

1.1 130E-18 

-1.1 130E-18 

Fujita[71] 




1.8300E-17 

0.00 


1.1 100E-18 

-1.1 100E-18 

Park[60] 




6.4200E-22 

0.49 


1.1200E-18 

-1.1200E-18 

Bird[ 1 5] 




2.6570E-13 

-0.98 


1.1 160E-18 

-1.1 160E-18 

Kang[64] 

39 

N + + N 2 - 

-> N 2 + + N 

1.6610E-18 

0.50 

1.0 

1.6840E-19 

-1.6840E-18 

Park[63] 




1.6700E-17 

-0.18 


1.6700E-19 

-1.6700E-19 

Bird[ 1 5] 




2.3250E-26 

2.10 


1.6700E-19 

-1.6700E-19 

Phelps[73] 




3.3540E-19 

0.81 


1.7950E-19 

-1.7950E-19 

Kang[64] 

40 

o + + n 2 - 

-> n 2 + + 0 

1.511 OE- 1 8 

0.36 

1.0 

3.1480E-19 

-3.1480E-19 

Park[74] 




1.4940E-18 

0.36 


3.1480E-19 

-3.1480E-19 

Losev[70] 




1.0550E-16 

-0.21 


3.0600E-19 

-3.0600E-19 

Bird[ 1 5] 




5.6460E-1 1 

-2.00 


3.1750E-19 

-3.1750E-19 

Kang[64] 

41 

0 + + NO - 

— > N + + 0 2 

2.3240E-25 

1.90 

1.0 

3.6730E-19 

-3.6730E-19 

Park[74] 

42 

N0 + + N - 

— > n 2 + + 0 

1.1960E-16 

0.00 

1.0 

4.9010E-19 

-4.1090E-19 

Park[74] 




2.8300E-17 

0.40 


4.9010E-19 

-4.9010E-19 

Bird[ 1 5] 

43 

N0 + + N - 

— > 0 f + n 2 

5.6460E-17 

-1.08 

1.0 

1.7670E-19 

-1.7670E-19 

Park[74] 

44 

N0 + + 0 - 

— > N + + 0 2 

1.6610E-18 

0.50 

1.0 

1.0660E-18 

-1.0660E-18 

Park[74] 




2.2250E-17 

0.31 


1.0670E-18 

-1.0670E-18 

Kang[64] 

45 

N0 + + 0 - 

— > 0 2 + + N 

1.1960E-17 

0.29 

1.0 

6.7100E-19 

-6.7100E-19 

Fujita[71] 

46 

NO 1 + 0 2 - 

-> 0 2 + + NO 

3.9850E-17 

0.41 

1.5 

5.4010E-19 

-5.4010E-19 

Park[74] 




3.9850E-17 

0.00 


4.5010E-19 

-4.5010E-19 

Losev[70] 




1.7200E-14 

-0.17 


4.4730E-19 

-4.4730E-19 

Bird[ 1 5] 




2.9890E-15 

0.17 


4.5560E-19 

-4.5560E-19 

Kang[64] 

47 

0 2 + + N - 

-> N + + 0 2 

1.4450E-16 

0.14 

1.0 

3.9490E-19 

-3.9490E-19 

Park[74] 


0/1 



Number 

Reaction 

a 

b 

C 

E a { J) 

E rx (J) 

Reference 

48 

o 2 + + N 2 - 

-*• n 2 + + 0 2 

1.6440E-17 

0.00 

1.5 

5.6190E-19 

-5.6190E-19 

Park[74] 

49 

0 2 + + 0 - 

■> 0 + + 0 2 

6.6420E-18 

-0.09 

1.0 

2.4850E-19 

-2.4850E-19 

Park[75] 




1.8900E-16 

-0.52 


2.5900E-19 

-2.5900E-19 

B ird [15] 




4.8490E-12 

-1.11 


3.8650E-19 

-3.8650E-19 

Kang[64] 


Table A. 15: Exothermic chemical reaction rate coefficients (m 3 /molecule/s) used in 
the TCE model. k r = al h cxp(-EJkT) 


Number 

Reaction 

a 

b 

C 

E a ( J) 

E rx (J) 

Reference 

50 

O + O + T 

->o 2 + t 

6.3050E-44 

-0.50 

0.0 

0.0000E+01 

8.1995E-19 

Farbar[6] 




8.3000E-45 

-0.50 


0.0000E+01 

8.1995E-19 

Gupta[76] 

51 

N + N + T 

-*n 2 + t 

5.6910E-40 

-1.50 

0.0 

0.0000E+01 

1.5630E-18 

Farbar[6] 




3.0060E-44 

-0.50 


0.0000E+01 

1.5630E-18 

Gupta[76] 

52 

N+O+T- 

-»• NO + T 

1.5830E-43 

-0.50 

0.0 

0.0000E+01 

1.0420E-18 

Farbar[6] 




2.7850E-40 

-1.50 


0.0000E+01 

1.0420E-18 

Gupta[76] 

53 

NO + N- 

■> n 2 + 0 

2.4900E-17 

0.10 

0.0 

0.0000E+01 

5.1750E-19 

Carlson[78] 




2.0200E-17 

0.10 


0.0000E+01 

5.9280E-19 

Bird[ 1 5] 




4.0600E-12 

-1.36 


0.0000E+01 

5.1750E-19 

Boyd[67] 




2.4900E-17 

0.00 


0.0000E+01 

5.1750E-19 

Gimelshein[65] 

54 

0 2 +N — > 

NO + O 

5.2000E-22 

1.29 

0.0 

4.9680E-20 

2.2380E-19 

Bird[ 15] 




4.1300E-21 

1.18 


0.0000E+01 

5.5300E-20 

Ozawa[9] 




4.6000E-15 

-0.55 


0.0000E+01 

2.2380E-20 

Boyd[67] 




1.5980E-18 

0.50 


4.9680E-20 

2.7190E-19 

Gimelshein[65] 

55 

N 2 + + e' - 

-> N +N 

4.4830E-12 

-0.50 

0.0 

0.0000E+01 

9.3400E-19 

Park[60] 




8.8800E-10 

-1.23 


0.0000E+01 

9.3400E-19 

Bird[ 1 5] 




1.0000E-08 

-1.43 


0.0000E+01 

9.3400E-19 

Ozawa[7] 




7.2740E-12 

-0.65 


0.0000E+01 

9.3400E-19 

Boyd[77] 




1.5700E-17 

0.85 


0.0000E+01 

9.3400E-19 

Boyd[67] 




2.4910E-08 

-1.50 


0.0000E+01 

9.3400E-19 

Kang[64] 

56 

NO + + e' ■ 

— > N + O 

1.4940E-11 

-0.65 

0.0 

0.0000E+01 

4.4220E-19 

Park[60] 




4.0300E-09 

-1.63 


0.0000E+01 

4.4220E-19 

Bird[ 1 5] 




2.0000E-06 

-2.05 


0.0000E+01 

4.4220E-19 

Ozawa[7] 




1.3210E-09 

-1.19 


0.0000E+01 

4.4220E-19 

Boyd[77] 




2.2000E-13 

-0.19 


0.0000E+01 

4.4220E-19 

Boyd[67] 




1.3280E-08 

-1.50 


0.0000E+01 

4.4220E-19 

Kang[64] 




2.9890E-1 1 

-1.00 


0.0000E+01 

4.4730E-19 

Gupta[72] 

57 

0 2 + + e' - 

-> O + O 

2.4900E-12 

-0.50 

0.0 

0.0000E+01 

1.1200E-18 

Park[60] 




3.8300E-09 

-1.51 


0.0000E+01 

1.1200E-18 

Bird[ 15] 




6.0000E-08 

-1.60 


0.0000E+01 

1.1200E-18 

Ozawa[7] 




1.4530E-04 

-2.41 


0.0000E+01 

1.1200E-18 

Boyd[77] 




9.2200E-15 

0.29 


0.0000E+01 

1.1200E-18 

Boyd[67] 




1.3280E-08 

-1.50 


0.0000E+01 

1.1 160E-18 

Gupta[72] 

58 

N 2 + + N — 

► N + + N2 

2.3700E-18 

-0.52 

1.0 

0.0000E+01 

1.6700E-19 

Bird[ 1 5] 




2.3430E-14 

-0.61 


0.0000E+01 

1.6700E-19 

Boyd[77] 




4.6500E-19 

0.50 


0.0000E+01 

1.6700E-19 

Kang[64] 

59 

N 2 " 1 " + 0 — 

> o f + n 2 

1.7700E-17 

-0.21 

1.0 

0.0000E+01 

3.1480E-19 

Bird[ 15] 




1.9780E-18 

0.11 


0.0000E+01 

3.1480E-19 

Boyd[77] 




4.1 180E-1 1 

-2.20 


0.0000E+01 

3.1480E-19 

Kang[64] 

60 

N + + 0 2 — > 

■ 0 + + NO 

2.4430E-26 

2.10 

1.0 

0.0000E+01 

2.1120E-19 

Boyd[77] 

61 

N 2 + + 0 — > 

• NO + + N 

4.1000E-18 

0.40 

1.0 

0.0000E+01 

4.9010E-19 

Bird[ 1 5] 











Number Reaction 

a 

b 

N + O + NO + — > NO + NO + 

0.150 

-1.90 

N + O + e' — > NO + e" 

0.000 

0.00 


Table A. 17: Coefficients required for the collision volume in Q-K exchange reactions 

(E a = a(T/T re j) h \E exch \; E a = E a - E exch if E exch < 0.0). 


Number Reaction 

34 N 2 + O -> NO + N 

35 NO + O — > 0 2 + N 

39 N + + N 2 — > N 2 + + N 

40 0 1 + N 2 — > N 2 + + O 

41 O +NO — > N + + 0 2 

42 NO + + N — > N 2 + + O 

43 NO + + N — >■ 0 + + N 2 

44 NO +0 — * N 4 + 0 2 

45 NO + + O -» 0 2 + + N 

46 NO + 0 2 — > OT + NO 

47 0 2 + + N — > N + + 0 2 

48 0 2 + + N 2 — » N 2 + + 0 2 

49 0 2 + + 0-^0 + + 0 2 

53 NO + N ->• N 2 + O 

54 0 2 + N — >NO + 0 

58 N 2 + + N — > N + + N 2 

59 N 2 + + O -> 0 + + N 2 

60 N + + 0 2 — > 0 1 + NO 

61 N 2 + + O — > NO + + N 

62 0 + + N 2 -> NO + + N 

63 N + + 0 2 — > NO 1 + O 

64 0 2 + + N — > NO 1 + O 

65 0 2 + + NO — > NO + + 0 2 

66 N 1 + 0 2 — > 0 2 + + N 

67 N 2 + + 0 2 — > 0 2 + + N 2 

68 0 + + 0 2 -^ 0 2 + +0 


a 

b 

Eexch (J) 

0.150 

0.00 

-5.2465E-19 

0.100 

0.68 

-2.6785E-19 

0.400 

0.10 

-1.6844E-19 

0.500 

0.20 

-3.1479E-19 

0.100 

0.50 

-3.6725E-19 

0.050 

0.00 

-4.9013E-19 

0.050 

0.70 

-1.7672E-19 

0.000 

0.00 

-1.0659E-19 

0.000 

0.00 

-6.7100E-19 

0.000 

0.00 

-4.5009E-19 

0.020 

0.10 

-3.9487E-19 

0.050 

0.55 

-5.6192E-19 

0.700 

0.15 

-2.4852E-19 

0.025 

0.90 

5.2465E-19 

0.100 

0.48 

2.6785E-19 

0.400 

0.25 

1.6844E-19 

0.600 

0.35 

3.1479E-19 

0.150 

0.40 

3.6725E-19 

0.020 

0.80 

4.9013E-19 

0.000 

0.00 

1.7672E-19 

0.000 

0.00 

1.0659E-18 

0.000 

0.00 

6.7100E-19 

0.070 

0.50 

5.4009E-19 

0.000 

0.00 

3.9487E-19 

0.200 

0.30 

5.6192E-19 

0.500 

0.10 

2.4852E-19 







Appendix B: Quantum-Kinetic Electronic Energy Level Transition Rates. 
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Figure B.l: Transitions for O2 + O2. 


Figure B. 4 : Transitions for O2 + N. 
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Figure B. 2 : Transitions for CE + N2. 


Figure B. 5 : Transitions for O2 + NO. 
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Figure B. 3 : Transitions for O2 + O. 
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Figure B. 6 : Transitions for N2 + O2. 
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Figure B.7: Transitions for N 2 + N 2 . 


Figure B. 10: Transitions for N 2 + NO. 
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Figure B.8: Transitions for N 2 + O. 
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Figure B. 1 1 : Transition from level 1-2 
for N 2 + e. 
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Figure B.9: Transitions for N 2 + N. 
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Figure B.12: Transition from level 2-3 
for N 2 + e. 
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Figure B.13: Transition from level 1-2 
for O + e. 
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Figure B.16: Transition from level 1-5 
for O + e. 
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Figure B.14: Transition from level 1-3 
for O + e. 
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Figure B.17: Transition from level 2-3 
for O + e. 
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Figure B.15: Transition from level 1-4 
for O + e. 
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Figure B.18: Transition from level 2-5 
for O + e. 
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Figure B.19: Transition from level 1-2 
for N + e. 



Figure B.22: Transition from level 1-5 
for N + e. 



Temperature (K) 

Figure B.20: Transition from level 1-3 
for N + e. 



Figure B.23: Transition from level 2-3 
for N + e. 



Figure B.21: Transition from level 1-4 
for N + e. 
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Figure B.24: Transition from level 2-4 
for N + e. 
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Figure B.25: Transition from level 2-5 
for N + e. 
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Figure B.28: Transitions for NO + O 2 . 
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Figure B.26: Transition from level 3-4 
for N + e. 
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Figure B.29: Transitions for NO + N 2 . 
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Figure B.27: Transition from level 3-5 
for N + e. 
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Figure B.30: Transitions for NO + O. 
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Figure B.3 1 : Transitions for NO + N. 
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Figure B.34: Transition from level 1-4 
for NO + e. 
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Figure B.32: Transitions for NO + NO. 
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Figure B.3 5: Transition from level 3-4 
for NO + e. 
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Figure B.33: Transition from level 1-3 
for NO + e. 
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Figure B.36: Transition from level 1-2 
for N 2 + + e. 
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Figure B.37: Transition from level 1-3 
for N 2 + + e. 
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Figure B.38: Transition from level 2-3 
for N 2 + + e. 


Table B. 1 : Sampled temperatures and their departure from the set equilibrium 
temperature. 


^equilibrium (K-) 

10000.00 

% Difference 

20000.00 

% Difference 

30000.00 

% Difference 

T,r(K) 

9949.81 

-0.502% 

19918.00 

-0.410% 

29901.50 

-0.328% 

T electron (K-) 

9953.79 

-0.462% 

19927.70 

-0.361% 

29942.00 

-0.193% 

Tel, 02 (K) 

10000.40 

0.004% 

20018.80 

0.094% 

30301.90 

1.006% 

Tel,N2 (K) 

10011.80 

0.118% 

20017.00 

0.085% 

30060.00 

0.200% 

T e i,o(K) 

10019.80 

0.198% 

19992.40 

-0.038% 

29967.50 

-0.108% 

Tci,n(K) 

10009.50 

0.095% 

19996.20 

-0.019% 

30040.60 

0.135% 

Tel, NO (K) 

10011.20 

0.112% 

20037.70 

0.189% 

30119.40 

0.398% 

Tel, 02+ (K) 

10012.50 

0.125% 

20046.50 

0.233% 

30140.30 

0.468% 

Tel, N2+ (K) 

10015.70 

0.157% 

20031.50 

0.158% 

30106.20 

0.354% 

Tel, 0+ (K) 

10011.30 

0.113% 

19993.80 

-0.031% 

29923.40 

-0.255% 

Tel,N+(K) 

10022.70 

0.227% 

19989.50 

-0.053% 

29896.20 

-0.346% 

T e i, no+ (K) 

10011.20 

0.112% 

20014.50 

0.073% 

30062.80 

0.209% 

Table B.2: 

Sampled and Boltzmann distributions for CF at 10,000 K. 



Level 

Sampled 

Boltzmann 

% 

Difference 

1 

7.7708E-01 

7.7719E-01 

-0.014% 

2 

1.6583E-01 

1.6585E-01 

-0.010% 

3 

3.8865E-02 

3.8807E-02 

0.149% 

4 

2.2448E-03 

2.2280E-03 

0.754% 

5 

1.0584E-02 

1.0566E-02 

0.175% 

6 

4.8032E-03 

4.7751E-03 

0.588% 

7 

5.9326E-04 

5.8615E-04 

1.212% 


1 0/1 


Table B.3: Sampled and Boltzmann distributions for O 2 at 20,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

5.5362E-01 

5.5315E-01 

0.086% 

2 

2.0892E-01 

2.0863E-01 

0.138% 

3 

7.1264E-02 

7.1363E-02 

-0.138% 

4 

1.7149E-02 

1.7099E-02 

0.292% 

5 

9.0780E-02 

9.1209E-02 

-0.470% 

6 

4.3126E-02 

4.3358E-02 

-0.535% 

7 

1.5139E-02 

1.5191E-02 

-0.342% 


Table B.4: Sampled and Boltzmann distributions for O 2 at 30,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

4.2546E-01 

4.2474E-01 

0.170% 

2 

1.9476E-01 

1.9370E-01 

0.550% 

3 

7.5215E-02 

7.5191E-02 

0.032% 

4 

2.9078E-02 

2.9007E-02 

0.245% 

5 

1.6010E-01 

1.6091E-01 

-0.506% 

6 

7.7062E-02 

7.7793E-02 

-0.939% 

7 

3.8326E-02 

3.8662E-02 

-0.869% 


Table B.5: Sampled and Boltzmann distributions for N 2 at 10,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

9.9500E-01 

9.9505E-01 

-0.005% 

2 

2.1963E-03 

2.1776E-03 

0.859% 

3 

1.1353E-03 

1.1238E-03 

1.026% 

4 

1.1045E-03 

1.0937E-03 

0.992% 

5 

2.1789E-04 

2.1573E-04 

1.000% 

6 

5.5589E-05 

5.4867E-05 

1.316% 

7 

9.4670E-05 

9.3263E-05 

1.509% 


Table B.6: Sampled and Boltzmann distributions for N 2 at 20,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

7.3201E-01 

7.3248E-01 

-0.064% 

2 

5.9496E-02 

5.9350E-02 

0.246% 

3 

6.0508E-02 

6.0296E-02 

0.352% 

4 

5.9591E-02 

5.9482E-02 

0.183% 

5 

1.8720E-02 

1.8681E-02 

0.211% 

6 

5.4476E-03 

5.4391E-03 

0.156% 

7 

1.0032E-02 

1.0029E-02 

0.034% 


1 















Table B.7: Sampled and Boltzmann distribution for N 2 at 30,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

3.6019E-01 

3.6023E-01 

-0.011% 

2 

9.7743E-02 

9.7283E-02 

0.473% 

3 

1.2460E-01 

1.2387E-01 

0.592% 

4 

1.2323E-01 

1.2275E-01 

0.391% 

5 

4.5125E-02 

4.5014E-02 

0.246% 

6 

1.371 IE-02 

1.3711E-02 

0.001% 

7 

2.5906E-02 

2.5975E-02 

-0.265% 


Table B.8: Sampled and Boltzmann distributions for O at 10,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

9.4529E-01 

9.4553E-01 

-0.025% 

2 

5.3797E-02 

5.3567E-02 

0.430% 

3 

8.1975E-04 

8.1264E-04 

0.875% 

4 

1.3074E-05 

1.2912E-05 

1.258% 

5 

5.1522E-06 

5.01 19E-06 

2.799% 

6 

6.0933E-06 

6.0882E-06 

0.084% 

7 

2.7848E-06 

2.7387E-06 

1.682% 


Table B.9: Sampled and Boltzmann distributions for O at 20,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

7.4519E-01 

7.4541E-01 

-0.029% 

2 

1.3215E-01 

1.3224E-01 

-0.070% 

3 

7.2810E-03 

7.2842E-03 

-0.044% 

4 

2.0615E-03 

2.053 IE-03 

0.409% 

5 

9.9362E-04 

9.9083E-04 

0.282% 

6 

2.4553E-03 

2.4419E-03 

0.550% 

7 

1.2762E-03 

1.2686E-03 

0.597% 


Table B.10: Sampled and Boltzmann distributions for O at 30,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

3.0872E-01 

3.0795E-01 

0.250% 

2 

8.0065E-02 

7.9930E-02 

0.170% 

3 

6.7504E-03 

6.7670E-03 

-0.245% 

4 

4.9942E-03 

4.9743E-03 

0.399% 

5 

2.5918E-03 

2.5813E-03 

0.406% 

6 

8.1065E-03 

8.0535E-03 

0.658% 

7 

4.4226E-03 

4.3898E-03 

0.748% 
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Table B. 1 1 : Sampled and Boltzmann distributions for N at 10,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

8.4636E-01 

8.4671E-01 

-0.042% 

2 

1.3341E-01 

1.331 IE-01 

0.225% 

3 

2.0090E-02 

2.0036E-02 

0.271% 

4 

1.5808E-05 

1.5762E-05 

0.292% 

5 

5.2667E-06 

5.2246E-06 

0.807% 

6 

7.9995E-06 

7.9046E-06 

1.201% 

7 

6.2563E-07 

6.0151E-07 

4.011% 


Table B.12: Sampled and Boltzmann distributions for N at 20,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

4.3952E-01 

4.3890E-01 

0.141% 

2 

2.7547E-01 

2.7515E-01 

0.116% 

3 

8.2188E-02 

8.2688E-02 

-0.605% 

4 

3.2850E-03 

3.2799E-03 

0.155% 

5 

1.3379E-03 

1.3353E-03 

0.198% 

6 

2.3279E-03 

2.3227E-03 

0.223% 

7 

2.6346E-04 

2.6158E-04 

0.719% 


Table B.13: Sampled and Boltzmann distributions for N at 30,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

1.0682E-01 

1.0620E-01 

0.580% 

2 

1.0632E-01 

1.0558E-01 

0.697% 

3 

3.9934E-02 

3.9954E-02 

-0.051% 

4 

5.8929E-03 

5.8551E-03 

0.646% 

5 

2.5708E-03 

2.5524E-03 

0.722% 

6 

4.6881E-03 

4.6517E-03 

0.783% 

7 

6.0134E-04 

5.9704E-04 

0.720% 


Table B.14: Sampled and Boltzmann distributions for NO at 10,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

9.8591E-01 

9.8601E-01 

-0.010% 

2 

8.3771E-03 

8.3267E-03 

0.605% 

3 

8.9642E-04 

8.8993E-04 

0.730% 

4 

1.3534E-03 

1.3421E-03 

0.843% 

5 

9.9780E-04 

9.8763E-04 

1.030% 

6 

8.6417E-04 

8.5528E-04 

1.040% 

7 

5.5615E-04 

5.5025E-04 

1.073% 
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Table B.15: Sampled and Boltzmann distributions for NO at 20,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

7.3783E-01 

7.3899E-01 

-0.157% 

2 

9.6384E-02 

9.6040E-02 

0.359% 

3 

1.5806E-02 

1.5699E-02 

0.684% 

4 

2.7450E-02 

2.7264E-02 

0.683% 

5 

2.3546E-02 

2.3388E-02 

0.675% 

6 

2.1886E-02 

2.1765E-02 

0.557% 

7 

1.7534E-02 

1.7457E-02 

0.439% 


Table B.16: Sampled and Boltzmann distributions for NO at 30,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

4.7041E-01 

4.7202E-01 

-0.342% 

2 

1.5317E-01 

1.5259E-01 

0.383% 

3 

2.9009E-02 

2.8737E-02 

0.948% 

4 

5.2855E-02 

5.2311E-02 

1.040% 

5 

4.7618E-02 

4.7228E-02 

0.826% 

6 

4.5222E-02 

4.5016E-02 

0.457% 

7 

3.8913E-02 

3.8862E-02 

0.132% 


Table B.17: Sampled and Boltzmann distributions for 02 + at 10,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

9.7135E-01 

9.7155E-01 

-0.020% 

2 

1.7156E-02 

1.7058E-02 

0.577% 

3 

2.8567E-03 

2.8334E-03 

0.821% 

4 

4.2667E-03 

4.233 IE-03 

0.795% 

5 

1.9302E-03 

1.9136E-03 

0.867% 

6 

5.8441E-04 

5.7826E-04 

1.063% 

7 

7.9150E-04 

7.8423E-04 

0.927% 


Table B.18: Sampled and Boltzmann distributions for 0? + at 20,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

6.6647E-01 

6.6768E-01 

-0.181% 

2 

1.2558E-01 

1.2512E-01 

0.372% 

3 

3.6258E-02 

3.6057E-02 

0.557% 

4 

5.4176E-02 

5.3977E-02 

0.369% 

5 

2.9739E-02 

2.9632E-02 

0.361% 

6 

1.1570E-02 

1.1518E-02 

0.451% 

7 

1.9025E-02 

1.8970E-02 

0.293% 
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Table B.19: Sampled and Boltzmann distributions for CL + at 30,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

4.1 158E-01 

4.1287E-01 

-0.313% 

2 

1.7107E-01 

1.7034E-01 

0.427% 

3 

5.9420E-02 

5.8989E-02 

0.731% 

4 

8.8882E-02 

8.8364E-02 

0.586% 

5 

5.1954E-02 

5.1754E-02 

0.386% 

6 

2.1926E-02 

2.1879E-02 

0.217% 

7 

3.8413E-02 

3.8443E-02 

-0.077% 


Table B.20: Sampled and Boltzmann distributions for N? + at 10,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

6.1227E-01 

6.1286E-01 

-0.096% 

2 

3.2813E-01 

3.2776E-01 

0.113% 

3 

1.5790E-02 

1.5718E-02 

0.457% 

4 

3.1316E-02 

3.1262E-02 

0.173% 

5 

6.2032E-03 

6.1665E-03 

0.596% 

6 

3.2954E-03 

3.2741E-03 

0.649% 

7 

6.1995E-04 

6.1403E-04 

0.964% 


Table B.2 1 : Sampled and Boltzmann distributions for N 2 + at 20,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

3.1457E-01 

3.1438E-01 

0.059% 

2 

3.2567E-01 

3.2514E-01 

0.163% 

3 

5.0346E-02 

5.0348E-02 

-0.003% 

4 

1.0017E-01 

1.0042E-01 

-0.245% 

5 

6.3131E-02 

6.3071E-02 

0.096% 

6 

4.5969E-02 

4.5958E-02 

0.025% 

7 

1.4075E-02 

1.4073E-02 

0.014% 


Table B.22: Sampled and Boltzmann distributions for N? + at 30,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

1.8259E-01 

1.8224E-01 

0.195% 

2 

2.3563E-01 

2.3481E-01 

0.349% 

3 

5.3865E-02 

5.3742E-02 

0.229% 

4 

1.0748E-01 

1.0728E-01 

0.183% 

5 

9.9579E-02 

9.9136E-02 

0.447% 

6 

8.0402E-02 

8.0276E-02 

0.157% 

7 

2.8918E-02 

2.8946E-02 

-0.098% 
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Table B.23: Sampled and Boltzmann distributions for 0 + at 10,000 K. 


Level Sampled Boltzmann % 

Difference 

~ 9.4566E-01 9.4590E-01 -0.025% 

2 5.01 10E-02 4.9904E-02 0.414% 

3 4.2276E-03 4.2007E-03 0.641% 


Table B.24: Sampled and Boltzmann distributions for 0 + at 20,000 K. 


Level Sampled Boltzmann % 

Difference 

1 6.9240E-01 6.9214E-01 0.037% 

2 2.5 1 3 IE-01 2.5137E-01 -0.023% 

3 5.6287E-02 5.6491E-02 -0.361% 


Table B.25: Sampled and Boltzmann distributions for 0 + at 30,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

5.2579E-01 

5.2459E-01 

0.229% 

2 

3.6206E-01 

3.6242E-01 

-0.099% 

3 

1.121 6E-01 

1.1299E-01 

-0.735% 


Table B.26: Sampled and Boltzmann distributions for N + at 10,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

9.4103E-01 

9.4 1 3 IE-01 

-0.030% 

2 

5.8007E-02 

5.7735E-02 

0.471% 

3 

9.5709E-04 

9.4864E-04 

0.891% 

4 

2.7276E-06 

2.7013E-06 

0.974% 

5 

1.3391E-07 

1.4098E-07 

-5.016% 

6 

0.0000E+01 

5.1 142E-10 

-100.000% 

7 

9.9063E-10 

6.0919E-10 

62.615% 


Table B.27: Sampled and Boltzmann distributions for N + at 20,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

8.3493E-01 

8.3487E-01 

0.008% 

2 

1.5404E-01 

1.541 IE-01 

-0.046% 

3 

8.8414E-03 

8.8344E-03 

0.079% 

4 

1.8248E-03 

1.8258E-03 

-0.056% 

5 

3.2492E-04 

3.2310E-04 

0.565% 

6 

1.4622E-05 

1.4505E-05 

0.810% 

7 

2.4824E-05 

2.4524E-05 

1.222% 
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Table B.28: Sampled and Boltzmann distributions for N + at 30,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

7.6006E-01 

7.5955E-01 

0.067% 

2 

2.0205E-01 

2.0243E-01 

-0.188% 

3 

1.7556E-02 

1.7600E-02 

-0.251% 

4 

1.5094E-02 

1.5173E-02 

-0.519% 

5 

4.0278E-03 

4.0337E-03 

-0.145% 

6 

4.1872E-04 

4.1885E-04 

-0.031% 

7 

7.9535E-04 

7.9589E-04 

-0.068% 


Table B.29: Sampled and Boltzmann distributions for NO + at 10,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

9.9599E-01 

9.9603E-01 

-0.004% 

2 

1.6518E-03 

1.6379E-03 

0.849% 

3 

1.1971E-03 

1.1854E-03 

0.988% 

4 

8.2190E-04 

8.1217E-04 

1.198% 

5 

1.7649E-04 

1.7472E-04 

1.015% 

6 

4.5549E-05 

4.4971E-05 

1.285% 

7 

6.9164E-05 

6.8332E-05 

1.218% 


Table B.30: Sampled and Boltzmann distributions for NO + at 20,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

7.8349E-01 

7.8358E-01 

-0.012% 

2 

5.5182E-02 

5.5037E-02 

0.263% 

3 

6.6323E-02 

6.6215E-02 

0.163% 

4 

5.4691E-02 

5.4809E-02 

-0.215% 

5 

1.7924E-02 

1.7975E-02 

-0.286% 

6 

5.2648E-03 

5.2652E-03 

-0.008% 

7 

9.1870E-03 

9.1786E-03 

0.091% 


Table B.3 1 : Sampled and Boltzmann distributions for NO + at 30,000 K. 


Level 

Sampled 

Boltzmann 

% 

Difference 

1 

4.5955E-01 

4.5873E-01 

0.178% 

2 

1.1342E-01 

1.1263E-01 

0.703% 

3 

1.6075E-01 

1.6052E-01 

0.145% 

4 

1.4048E-01 

1.4151E-01 

-0.728% 

5 

5.2954E-02 

5.3415E-02 

-0.862% 

6 

1.6242E-02 

1.6335E-02 

-0.568% 

7 

2.9677E-02 

2.9810E-02 

-0.446% 
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Appendix C: Quantum-Kinetic Reaction Rates. 
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Figure C. 1 : Q-K reaction rates for 
reaction 1. 
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Figure C.3: Q-K reaction rates for 
reaction 3. 
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Figure C.2: Q-K reaction rates for 
reaction 2. 



Figure C.4: Q-K reaction rates for 
reaction 4. 
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Figure C.5: Q-K reaction rates for 
reaction 5. 
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Figure C.8: Q-K reaction rates for 
reaction 8. 
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Figure C.6: Q-K reaction rates for 
reaction 6. 
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Figure C.9: Q-K reaction rates for 
reaction 9. 
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Figure C.7: Q-K reaction rates for 
reaction 7. 
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Figure C. 10: Q-K reaction rates for 
reaction 10. 


1 1 t. 


1 i i r 

N +0 > N + N + 0 



0 5xl0 3 lxio 4 2xl0 4 2xl0 4 3xl0 4 3xl0 4 4xl0 4 

Temperature (K) 

Figure C. 1 1 : Q-K reaction rates for 
reaction 1 1 . 
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Figure C. 14: Q-K reaction rates for 
reaction 14. 
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Figure C. 12: Q-K reaction rates for 
reaction 12. 
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Figure C. 15: Q-K reaction rates for 
reaction 15. 
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Figure C. 13: Q-K reaction rates for 
reaction 13. 
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Figure C. 16: Q-K reaction rates for 
reaction 16. 


1 1 A 



Figure C. 17: Q-K reaction rates for 
reaction 17. 
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Figure C. 1 8: Q-K reaction rates for 
reaction 18. 
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Figure C.20: Q-K reaction rates for 
reaction 20. 
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Figure C.2 1 : Q-K reaction rates for 
reaction 2 1 . 
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Figure C. 19: Q-K reaction rates for 
reaction 19. 
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Figure C.22: Q-K reaction rates for 
reaction 22. 
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Figure C.23: Q-K reaction rates for 
reaction 23. 
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Figure C.26: Q-K reaction rates for 
reaction 26. 
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Figure C.24: Q-K reaction rates for 
reaction 24. 
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Figure C.27: Q-K reaction rates for 
reaction 27. 
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Figure C.25: Q-K reaction rates for 
reaction 25. 
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Figure C.28: Q-K reaction rates for 
reaction 28. 
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Figure C.29: Q-K reaction rates for 
reaction 29. 



Figure C.32: Q-K reaction rates for 
reaction 32. 
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Figure C.30: Q-K reaction rates for 
reaction 30. 
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Figure C.33: Q-K reaction rates for 
reaction 33. 
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Figure C.3 1 : Q-K reaction rates for 
reaction 3 1 . 
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Figure C.34: Q-K reaction rates for 
reaction 34. 
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Figure C.35: Q-K reaction rates for 
reaction 35. 
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Figure C.38: Q-K reaction rates for 
reaction 38. 
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Figure C.36: Q-K reaction rates for 
reaction 36. 
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Figure C.39: Q-K reaction rates for 
reaction 39. 
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Figure C.37: Q-K reaction rates for 
reaction 37. 
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Figure C.40: Q-K reaction rates for 
reaction 40. 
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Figure C.4 1 : Q-K reaction rates for 
reaction 4 1 . 
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Figure C.44: Q-K reaction rates for 
reaction 44. 



Temperature (K) 

Figure C.42: Q-K reaction rates for 
reaction 42. 



Figure C.45: Q-K reaction rates for 
reaction 45. 



Figure C.43: Q-K reaction rates for 
reaction 43. 
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Figure C.46: Q-K reaction rates for 
reaction 46. 
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Figure C.47: Q-K reaction rates for 
reaction 47. 



E 

CS 

■4-J 

<fl 

£ 

O 

o 


<D 

re 


DL 


IQ ’ 44 


10' 45 


10- 46 


10' 47 

0 5000 10000 15000 20000 25000 30000 35000 




• 

Q-K Sampled 

0 

— Park[65]* 

B 

— Farbar[6] 

e — 

— Gupta[78] 



Temperature (K) 

Figure C.50: Q-K reaction rates for 
reaction 50. 
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Figure C.48: Q-K reaction rates for 
reaction 48. 
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Figure C.5 1 : Q-K reaction rates for 
reaction 5 1 . 



Figure C.49: Q-K reaction rates for 
reaction 49. 
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Figure C.52: Q-K reaction rates for 
reaction 52. 
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Figure C.53: Q-K reaction rates for 
reaction 53. 
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Figure C.56: Q-K reaction rates for 
reaction 56. 
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Figure C.54: Q-K reaction rates for 
reaction 54. 


-2 

a> 

3 

O 

o 

o 

E 


c 

to 

</> 

c 

o 

o 


to 

a: 



Temperature (K) 

Figure C.57: Q-K reaction rates for 
reaction 57. 
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Figure C.55: Q-K reaction rates for 
reaction 55. 
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Figure C.58: Q-K reaction rates for 
reaction 58. 
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Figure C.59: Q-K reaction rates for 
reaction 59. 
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Figure C.60: Q-K reaction rates for 
reaction 60. 



Figure C.6 1 : Q-K reaction rates for 
reaction 6 1 . 



Figure C.62: Q-K reaction rates for 
reaction 62. 



Figure C.63: Q-K reaction rates for 
reaction 63. 
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Figure C.64: Q-K reaction rates for 
reaction 64. 
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Figure C.65: Q-K reaction rates for 
reaction 65. 
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Figure C.68: Q-K reaction rates for 
reaction 68. 
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Figure C.66: Q-K reaction rates for 
reaction 66. 



Figure C.67: Q-K reaction rates for 
reaction 67. 
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